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The June 1992 to May 1993 grant NCC-2-677 provided for the 
continued demonstration of Computational Fluid Dynamics (CFD) as 
applied to the Stratospheric Observatory for Infrared Astronomy 
(SOFIA). While earlier grant years allowed validation of CFD through 
comparison against experiments, this year a new design proposal was 
evaluated. The new configuration would place the cavity aft of the 
wing, as opposed to the earlier baseline which was located 
immediately aft of the cockpit. This aft cavity placement allows for 
simplified structural and aircraft modification requirements, thus 
lowering the program cost of this national astronomy resource. 

The configuration, described in Appendices C and D, was 
computed at 41,000 foot cruise altitude conditions. The aft-ramp 
cavity treatment demonstrated similar broadband acoustic quieting 
behavior as seen in the earlier baseline designs. However, although 
the shear layer reattaches along the aft ramp, the low-momentum 
flow separates behind the cavity. This separation is of considerable 
concern due to its proximity to the tail control surfaces. The 747SP 
horizontal and vertical stabilizer geometry, only recently received 
from Boeing, were included in a preliminary computation. If 
separation about the full geometry is found to be significant, then a 
sizeable part of the time and expense of a wind tunnel test and 
model would have been saved. 

In addition to studying the effect of configuration changes on 
the flow, the effect on optical clarity was also investigated, as 
documented in Appendix D. The computed flow about two U.S. Army 
Airborne Optical Adjunct (AOA) configurations provided validation, 
the resultant methodology was then applied to the aft cavity SOFIA 
design. Tentative comparison of forward and aft cavity placement 
shows some degradation of clarity in the near-infrared for the aft 
location. 

The past year has seen the continuing design of SOFIA, and CFD 
has begun to provide information, such as telescope loadings and 
flow separation regions, which are difficult and expensive to 
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experimental means. The generality of the developed methodology 
has attracted industry interest for applications ranging from the 
stealthy internal carriage of stores to the design of airborne laser 
platforms. 

Numerical evaluation of the fluid and optical fields will provide 
a means, in addition to experimental results, to evaluate 
configurations in a timely and cost effective way. The results of this 
work will provide a framework for industrial design of 
configurations which experience related flow phenomena. 

In addition, this grant period saw the commencement of an 
investigation into the coupling of fluids, rigid-body dynamics, and 
active controls. Appendix E summarizes the demonstration of the 
coupled system for a two-dimensional application. The simulation 
was performed on a vector supercomputer, the Numerical 
Aerodynamic Simulator (NAS) Cray Y-MP, and was validated against 
linearized dynamics using analytic sensitivity coefficients. The 
fourth-order system for pitch attitude control used proportional and 
rate gyro feedback, and included a canard servo model. The gains 
were determined using conventional pole placement techniques, and 
were held fixed during the course of the nonlinear coupled 
simulations. 

The coupling of an active control system with nonlinear fluid 
and body dynamics will allow simulation of aircraft in critical 
regimes of the flight envelope. This computational prototyping offers 
reduced development costs while enhancing safety for the next 
generation of high performance aircraft. 
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APPENDIX A 



Abstract 


Numerical simulation of two classes of unsteady flows are obtained via the Navier- 
Stokes equations: a blast-wave/target interaction problem class and a transonic cavity 
flow problem class. The method developed for the viscous blast-wave/target interac- 
tion problem assumes a laminar, perfect gas implemented in a structured finite- volume 
framework. The approximately factored implicit scheme uses Newton subiterations to 
obtain the spatially and temporally second-order accurate time history of the interac- 
tion of blast-waves with stationary targets. The inviscid flux is evaluated using either 
of two upwind techniques, while the full viscous terms are computed by central differ- 
encing. Comparisons of unsteady numerical, analytical, and experimental results are 
made in tw r o- and three-dimensions for Couette flows, a starting shock-tunnel, and 
a shock-tube blockage study. The results show accurate wave speed resolution and 
nonoscillatory discontinuity capturing of the predominantly inviscid flows. Viscous 
effects were increasingly significant at large post-interaction times. 

While the blast-wave/target interaction problem benefits from high- resolution 
methods applied to the Euler terms, the transonic cavity flow problem requires the 
use of an efficient scheme implemented in a geometrically flexible overset mesh envi- 
ronment. Hence, the Reynolds averaged Navier-Stokes equations implemented in a 
diagonal form are applied to the cavity flow class of problems. Comparisons between 
numerical and experimental results are made in two-dimensions for free shear layers 
and both rectangular and quieted cavities, and in three-dimensions for Stratospheric 
Observatory For Infrared Astronomy (SOFIA) geometries. The acoustic behavior of 
the rectangular and three-dimensional cavity flows compare well with experiment in 
terms of frequency, magnitude, and quieting trends. However, there is a more rapid 
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decrease in computed acoustic energy with frequency than observed experimentally 
owing to numerical dissipation. In addition, optical phase distortion due to the time- 
varying density field is modelled using geometrical constructs. The computed optical 
distortion trends compare with the experimentally inferred result, but underpredicts 
the fluctuating phase difference magnitude. 
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Chapter 1 
Introduction 


Performance envelopes of vehicles and structures which interact with the atmosphere 
are often limited by unsteady aerodynamic effects. Specific examples which are ad- 
dressed here range from the blast- wave/target interaction problem, where peak over- 
pressures are many times quiescent conditions, to the seeing problem, where densitv 
fluctuations contribute to image degradation. Only recently have advances in com- 
puting power and numerical algorithms provided the potential, complementary to 
experimental studies, for the timely design of effective configurations. However, the 
use of unsteady computations in the design phase is presently at an immature stage 
of development. The objective of this effort is to demonstrate computational tech- 
nologies as applied to current topics of interest in the unsteady, compressible perfect 
gas regime. It is hoped that, through comparison to accepted experimental data, 
computational methods can make significant contributions to the design of systems 
which interact with unsteady flows. 

In these studies, two solution methods to the Navier-Stokes equations are pre- 
sented and applied to several test cases. The cases pertain to the blast- wave/ target 
interaction or the transonic cavity flow problem classes. The method used to model 
the strongly unsteady blast-wave flows concentrates on resolution of the complex 
physics through the use of characteristic-based schemes. In contrast, for the tran- 
sonic cavity flow problem, the combination of complex geometries and large problem 
size requires the use of an efficient integration scheme. Comparison of the numerical 
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and experimental results will provide a reference base of unsteady numerical results 
for use in configuration design. 

1.1.1 Blast- Wave Problem 

The study of the effects of blast-wave impingement upon vehicles and structures is 
of practical consideration in the determination of their survivability. The experimen- 
tal study of the blast- wave/target interaction problem requires the use of expensive 
above ground tests or facilities such as the U.S. Army Large Blast/Thermal Simulator 
(LB/TS) facility depicted in Fig. 1. Moreover, experiments can suffer from limited 



Figure 1: Proposed Large Blast/Thermal Simulator facility [1] 


phase durations, and deduction of the physics of the flowfield is difficult because 
of practicalities in data acquisition methods. Visualization of the propagating wave- 
fronts allows separation of pressure peaks due to wave reflection from shock tube walls 
from the pressure history. These pressure spikes can then, for example, be removed 
from structural frequency excitation analyses. Therefore, the information obtained 
via numerical simulation can be used for design from dynamic similitude conditions, 
and to augment data obtained in the test facilities. 

The simulation of the blast-wave problem has been studied in varying degrees 
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of physical complexity, from self-similar Euler to two-dimensional viscous flows. The 
simulation of the blast-wave problem using the Euler equations was studied by Kutler. 
Sakell, and Aiello [2] in 1975, where the intersection of a planar shock with a wedge 
in supersonic flight was modelled using the MacCormack scheme [3]. The self-similar 
nature of these types of flow with respect to time was used to obtain two- and later 
three-dimensional cone solutions [4]. In addition, Kutler and Shankar [5] used a shock- 
fitting procedure for the regular diffraction of a planar shock by a wedge, which is 
also self-similar. Although good results can be obtained using discontinuity fitting 
methods, coding complexities have generally hindered their application for complex 
situations. 

A general solution of the truly time-dependent inviscid interaction problem was 
modelled by Champney, Chaussee, and Kutler [6] in 1982. Shock diffraction over sev- 
eral simple two- and three-dimensional geometries were presented using the 
MacCormack and Beam- Warming schemes [7]. Mark and Kutler [8] performed a 
two-dimensional simulation of a shock passing over a simplified profile of a truck. 
However, inaccuracies due to discontinuity smearing and oscillations led to the de- 
velopment and application of explicit and implicit high-resolution schemes in the 
mid- 1980’s. Several two-dimensional problems using these high-resolution methods 
were presented by Yee [9] and Hisley and Molvik [10]. Recently, Lohner [11] has ob- 
tained solutions for complex geometries via an adaptive unstructured approach. The 
geometric flexibility of an unstructured grid method offers great potential, however 
the use of stretched tetrahedral grids required for efficient computation of viscous 
flows is a current topic of research. 

The solution of the Navier-Stokes equations in the high Reynolds number regime 
requires the use of fine grids to resolve thin viscous layers. The concomitant stiff- 
ness arising from the difference in length scales suggests the use of implicit schemes. 
Bennett, Abbett, and Wolf [12] applied a Beam- Warming scheme [7] with the Baldwin- 
Lomax [13] turbulence model to the problems of a developing boundary-layer behind 
a moving shock and shock diffraction over a cylinder. Molvik [14] used implicit 
high-resolution methods to obtain two-dimensional solutions for the unsteady devel- 
opment of a boundary-layer, the cylinder diffraction problem, and intersection of a 
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planar shock with a missile in supersonic flight. 

The purpose of this effort is to address the general three-dimensional, viscous 
blast-wave problem. The techniques developed here utilize total variation diminishing 
(TVD) upwind and upwind-biased schemes to resolve the discontinuous flow features 
without the oscillations prevalent in the more conventional central difference methods. 
Wave speeds are resolved adequately at large Courant numbers through the use of 
time conservative differencing and Newton subiterations. 


1.1.2 Cavity Flow Effort 

The Stratospheric Observatory For Infrared Astronomy (SOFIA) will be a three meter 
class Cassegrain telescope which utilizes a Boeing 747SP as an observation platform. 
An artist’s concept of the observatory, which is a follow-on to the Kuiper Airborne 
Observatory, is shown in Fig. 2. This airborne system, currently being studied by the 



Figure 2: Artist's concept of the SOFIA configuration 


United States NASA and Germany’s DARA, offers capabilities which augment land 
and space-based options in several ways. First, the mission flexibility of a long-range 
mobile platform lends astronomers freedom to investigate transient astronomical phe- 
nomena on a global basis. Second, atmospheric attenuation of some wavelengths 
of interest provide motivation for a platform which operates above the tropopause. 
Third, the cost of maintaining and upgrading observation technologies is lower than 
would be incurred with an orbiting configuration. 
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Nevertheless, the use of an aircraft-based observatory presents some challenges. 
The limited bandwidth of solid materials in the infrared frequency ranges of inter- 
est preclude their use as windows. The temperature of the window material would 
also contribute to background radiation levels. Therefore, the telescope cavity must 
remain open to the freestream. Empirical evidence has shown that violent shear 
layer oscillations with concomitantly dangerous levels of acoustic loading occur for 
untreated, or rectangularly shaped, open cavity configurations [15, 16]. Hence, there 
is a need to develop cavity flow control treatments to suppress the flow unsteadiness, 
both to reduce the risk of injury to the crew and to obtain high-quality seeing. To- 
wards these objectives, both experimental and computational fluid dynamics (EFD 
and CFD) analyses will be used in the design cycle. The purpose of this work is 
to develop and apply numerical tools for use in the design of the next generation 
airborne observatory. 

The driven cavity problem has been a subject of much research, both experimental 
[15-29] and numerical [30-40], owing to its wide practical applicability. The buffeting 
and sound production of bomb bays, slotted wind-tunnel walls, transition-delaying 
airfoil cavities, and deflected control surfaces are examples of the range of problems. 
The effort here is focused upon the transonic regime, and previous efforts in these 
flow speeds are reviewed by Komerath, Ajuha, and Chambers [41] and Rockwell and 
Naudascher [42]. Some of the research which is pertinent to the transonic aero-window 
problem is highlighted below. 

In 1955, Karamcheti [17] studied subsonic and low supersonic flow over rectangular 
cavities, in which the inverse relationship between cavity length, L, and dominant res- 
onant frequency, as well the acoustic intensity and radiation directivity w r as reported. 
In the same year, Roshko [18] documented skin friction and pressure distribution along 
the cavity walls. The three-dimensional low subsonic study of Maull and East [20] 
showed that spanwise cells can modestly perturb C p distributions from the idealized 
infinitely wide cavity. Rossiter [15], in 1964, related cavity resonance to the edge-tone 
phenomenon, and deduced a model applicable to the transonic regime of concern here. 
In this study, Rossiter also demonstrated that a cavity leading edge spoiler drasti- 
cally reduced acoustic levels. The work of Heller and Bliss [24] demonstrated that 
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the cavity feedback mechanism can also be suppressed by geometry modifications at 
the shear layer impingement region. These wind tunnel tests have provided valuable 
insight into the governing flow mechanisms, but dynamic similitude was typically not 
achieved. The dangerous behavior of cavity flows generally precludes the use of flight 

tests to verify scaling laws. Numerical simulation offers the potential for safely scaling 
sub-scale tests to flight conditions. 

Toward this objective, numerical efforts of the type which can model the viscous 
flowfield about a geometrically complicated structure was begun in 1979 by Hankey 
and Shang [30], Using MacCormack’s scheme, the dominant resonant mode of a rect- 
angular cavity was accurately predicted. Ora [37] used the same scheme to compute 
flow about quieted two-dimensional cavities. In 1987, Suhs [34] used a block im- 
plicit scheme to obtain the viscous flow about a parallelepiped cutout. The overset 
mesh method which was the precursor to that used herein was utilized. Dougherty 
et al. [39] have recently completed a detailed study of two-dimensional cavities using 
a high resolution scheme. They computed distinct spectral peaks which have been 
observed experimentally. 

The present effort builds on past numerical studies by validating the ability to 
predict free shear layers, and both untreated and treated two- and three-dimensional 
cavity configurations with an efficient scheme in an overset mesh framework. Compar- 
ison of computed flowfields to experimental and analytical results allows assessment 
of cavity load prediction capabilities. Prediction of the acoustic intensity levels and 
frequencies are of primary interest for safety reasons. However, estimation of optical 
distortion is required to determine mission effectiveness. 

1.1.3 Aero-Optics Work 

The study of the effect of a fluid field upon an optical field, dubbed aero-optics, has 
been extensive over the past four decades. Applications include imaging of re-entry 
vehicles or, as in this study, of astronomical bodies through the atmosphere. Many 
experimental and theoretical approaches to the optical distortion problem have been 
investigated, those of which are pertinent to this transonic aero- window problem are 
summarized here. The experimental efforts can be grouped into two categories: direct 
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measurement methods and techniques based on aerodynamically inferred quantities. 
Results obtained via the latter method are more prevalent because of practical diffi- 
culties in direct measurement techniques [43]. In fact, only aerodynamically inferred 
distortion levels will used for validation in the present work. 

Although early experimental and theoretical efforts assumed incoherent statistical 
turbulence [44, 45, 46], recent studies have begun to examine the effect of shear layer 
structures on electromagnetic field distortion. Using a passive scalar field from a 
direct numerical simulation, Truman and Lee [47] found an optimum viewing angle 
normal to the hairpin vortices in the homogeneous sheared fluid region. They also 
found analysis via non-refracting geometric optics to be equivalent to the parabolized 
Helmholtz representation of light. Although this class of studies provides excellent 
insight into the effects of small-scale structure on the electromagnetic field, it is clear 
that the expense of such methods precludes their near-term use for the problems 
under consideration here. 

The study of large scale structures in shear layers has been an active topic of 
research since they were observed by Brown and Roshko in 1974 [48]. Only recently 
has the effect of these structures on the optical field been studied. In 1990, Chew and 
Christiansen [49, 50] experimentally observed the effect of shear layer structures on 
beam propagation. Tsai and Christiansen [51] used an Euler simulation to determine 
the optical characteristics of a perturbed free shear layer. The use of a growing 
sinusoidal phase plate to represent the effect of vortical structures on an optical field 
was hypothesized. Wissler and Roshko [52] recently performed an experimental study 
of the motion of a thin light beam caused by passage through a shear layer. They 
postulated that spanwise steering asymptotes to a higher level than the streamwise 
component. 

The numerical modelling of the optical effect of a cavity-spanning shear layer was 
presented by Cassady et al. [53] in 1987. They found their two-dimensional solution 
to result in poor prediction of optical distortion. Farris and Clark [54, 55] used time- 
mean quantities and empirical evidence to ascertain the fluctuating density levels 
required for optical analysis. 
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The present effort attempts to determine what portion of the optical path distor- 
tion can be resolved using cell sizes required to obtain an accurate flowfield solution. 
Computed optical distortion levels are compared to flight and wind tunnel measure- 
ments for two- and three-dimensional quieted cavities, respectively. The following 
chapters address the methods used to predict the unsteady flows, the modelling of 
turbulence, and the optical distortion. 



Chapter 2 

Numerical Method 


2.1 The Governing Equations 


The Navier-Stokes equations may be expressed in integral conservation law form, 
coupled with the continuity and energy equations as 

I (1) 

where body forces have been neglected and the cell volumes are time invariant. Here 
V is the volume of an arbitrary fluid packet, F = Es r s‘ l + -fvsj + G/vsk is the flux 
tensor of second order, and ds is an outward directed normal of a differential surface 
area. The vectors may be written in Cartesian coordinates as 

Q = [p,pu,pv,pw,e] T 


Ens = 


pu 

pU 2 +p+ T xx 
PUV + T x y 
PUW + T xz 

(e + p)u + T XX U + T X yV + T X: W + q x 
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F ns 


pv 

pvu + T yx 

PV 2 +P+Ty y 
PVW + Ty z 

(e + p)v + Ty 2 U + TyyV + Ty ; W + Qy 


G\s = 


pw 

pwu + T zx 
pwv + T zy 
pw 2 + p + t 2Z 
(e + p)w + T ZX U + T Z yV + t Z2 w + q z J 


where each flux can be partitioned into inviscid and viscous portions. The density, 
pressure, and velocity components are respectively given by p,p,u,v, and w. The 
viscous stresses are composed of the terms: 


Trx 

T yy 


dv dw ' 
dy dz 


_ 9 du ( du 

2 ^dx A (ax 

n dv ( du dv dw \ 

- 2 %- x {T x + d-» + j;) 



The total energy per unit volume, e, is related to the internal energy per unit mass, 
e, by e = pe + pq 2 / 2. The perfect gas equation of state, p = pRT, completes the 
system. In addition, for thermally and calorically perfect gases, the internal energy 
per unit mass and the enthalpy per unit mass can be expressed solely as functions of 
temperature: 


dc = c v dT , dh = c p dT 
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Finally, for a calorically perfect gas the specific heats are constant, leaving e = c v T + 
const., and h = c p T + const., where the additive constants may be set to zero. The 
ratio of specific heats and specific gas constant are 


Cp i 


7 = — = - for air, R = c v — c v 
C v 5 

and the thermodynamic variables are related using 


/lT — 


( e +p) P , P, 2 , 2 , 

, e= - + -(u +v i + w i ) 

p 7-12 


where hj is the total enthalpy per unit mass. 

Fourier’s law for heat transfer by conduction is assumed; hence, the heat transfer 
can be expressed as 


q 


-kVT = - (q x i + q y j + ^k) 
f de. dc . de \ 

" K (al 1 + + dl k ) 


w’here k = k/c v = 7 p/Pr. The Prandtl number for air, which is a function only of 
the gas, relates the diffusion of momentum to the diffusion of heat, and is fixed at 
Pr = 0.72. 

The relationship between the first (p), second (A), and bulk (£) viscosity coeffi- 
cients is £ = i// + A. The bulk viscosity coefficient is set to zero in accordance with 
Stokes’ hypothesis, resulting in A = — |/r. This hypothesis is invoked here based on 
the assumption that the relative effects of the shearing stress is much larger than those 
caused by the dilational stress effects, not on the theory for monatomic gases [56]. 
By using the kinetic theory of gases, the physical phenomena of thermal conductivity 
and viscosity can be expressed in terms of the thermodynamic states. Viscosity is 
related to the thermodynamic state using Sutherland’s formula: 


P = 


CiTt 

r + c 2 


where C\ and C 2 are specific to the gas in question. 
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2.2 Turbulence Model 


Current limitations in computing power relegate most engineering computations to 
the use of grids which are too coarse to resolve all of the pertinent scales of motion. 
Reynolds averaging of the Navier-Stokes equations decomposes the flow into slowly 
and rapidly varying components [57], The slowly varying component is resolved from 
the spatial and time integration step sizes, while the rapidly fluctuating component 
is modelled. Although the blast-wave computations did not model turbulence, the 
turbulence model used in the cavity flow problems is outlined here for completeness. 

The effective or eddy viscosity due to additional turbulent mixing can be related 
to mean stresses using the Boussinesq approximation. The total effective viscosity is 
then given by ytotai y molecular y turbulent ~ y 4“ The Reynolds stress resulting 

from the Boussinesq assumption is 


-r-r (du, duj\ 2 ( dTT k \ 


where ( ) denotes a time mean. The eddy viscosity is given by p t oc piv where the 
length and velocity scales are given by £ and v. Alternatively, the expression for the 
eddy viscosity given by Prandtl is p t oc pP\u\, where the magnitude of vorticity is 


u = 


i 


^ du dv ^ 2 + / dv dw ^ 2 ^ ( dw du 


\dy dx t 


dy) 


dx di 


The local density is specified by Morkovin’s hypothesis, which states that compress- 
ibility does not affect the scales of the turbulent motion. 

The algebraic turbulence model of Baldwin and Lomax [13], as modified and 
implemented by Buning [58], is described below. This description is included to clearly 
show how the modified model constant used in the computed cavity shear layers is 
determined. The description assumes flow in the ( x,y ) plane with the freestream 
aligned with the x coordinate. 


Treatment of Wall Bounded Flows 

In Prandtl’s mixing length hypothesis the turbulent eddy size is limited by the prox- 
imity to the wall, giving I - ky , where the von Karman constant k = 0.4. The 
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addition of Van Driest wall damping results in 


t = ky(l - e - y+/4+ ) 

where .4 + = 26, y + = r “ - , and the wall shear stress is t w = (^f^) • 

The Baldwin-Lomax two-layer model uses the Prandtl-Van Driest model for the 
inner layer, and is given by 


{ pp M y < Vcrosaover Prandtl-Van Driest 

pC cpR EwakeE Kleb V ^ Vcrossove r Outer region 

where y is the distance from the wall, ycrossover is the location of the first intersection 
of inner and outer values of p t , the Clauser constant is K = 0.0168, C ^ — 1.6, and 


wake 


— min | 


VmaxF i 

CwkVi 


max 


max p 

‘mu 


where C w k = 1.0. The quantities F max and its location y max are found from 


F{y) = y|w| (l - e 


where the exponential term is dropped in wakes. The search for the Fmax term ends 
when F(y) drops below a specified percent of the first peak away from the wall, or at 
overset mesh boundaries. The Klebanoff intermittency function is specified according 
to 


F k'leb(y) = 


1 + 5.5 


CjOeby\ 

Vmax ) 


where CKieb = 0.3. The total velocity difference is given by 


ujif = \J (a 2 + v 2 )^ - v /(u 2 + u 2 ) r 


and the latter term is taken as zero at the wall. 


Treatment of Free Shear Flows 

The eddy viscosity in the shear layer was computed as outlined by Buning [58]. 
Development of the free shear layer model begins by using F(y) = y|u>|, as suggested 
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by Baldwin and Lomax [13] for wake regions. This results in 

,2 


wake 


— C w lc 


ymax^dif 


= ^frr^) 

VMmar/ 

oc C 2 \uj\ 


M. 


where specification of C w k is discussed below and the velocity difference is modified 
to be half the total velocity difference between the streams in the specified shear layer 
region 

Udif = v /(« 2 + V 2 ) - t /(||2 + v 2 ) 

v 'max V V 

Finally, the Klebanoff intermittency function is modified to 


Fh'iebiy ) = 


i + 5.5 f CKi *\y ~ yM».» 

V Vw 


61 


-1 


where the shear layer width is given by y w = u dif / \u\ max . The free shear layer model 
is now given by 

Hi — pKCcpCukr-j^— ( 2 ) 

Mmar V ’ 

after dropping the intermittency function for the analysis below. 

The magnitude of the eddy viscosity in the free shear layer model can be altered 
by specification of C w k. The remainder of this section shows how C w h is chosen based 
on empirical and analytic information. 

Gortler’s shear layer solution is given by 


u = 


U\ A U 2 


1 + 


«2 — Wj 

U2 + 1*1 


er/(r?) , erf(rj) = ~^= f e * drj , 77 

Y 7T JO 


<?y 

X 


where u\ and u 2 are the velocities of the slow and fast streams and 77 is the similarity 
coordinate (See Fig. 19). The spreading parameter <r is inversely related to the 
spreading rate, db/dx, where b is a measure of the shear layer width. The value of 
the spreading parameter when the velocity of one of the streams is zero is < 7 q. 
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Gortlers solution can be used to determine the maximum vorticity magnitude in 
a free shear layer as follows: 

du dr} 
dr} dy 

: ie~ n 2 , A u = («2 — u i) 

Mmax = —^Au 
XpK 

Now, using Prandtl’s mixing length assumption and scaling laws for jet bound- 
aries, eddy viscosity can also be expressed as 



Pt 


OC pl 2 \ui\ 

2 Au 

pl T 

= KopbAu 


( 3 ) 


where K 0 = ^ t- 

Setting Eqs. 2 and 3 equal results in 

C wk = (aoKC^r 1 

and only cto remains to be specified. Estimates of <Jo from empirical evidence is quite 
variable, ranging from 9.0 to 13.5 primarily dependent upon whether the upstream 
boundary layer is turbulent or laminar [48, 59, 60]. For this series of cavity flow 
efforts (Tq was set to 11.0, which appears to be the result from the highest quality 
experiments, resulting in a value of C wk = 1.91. Previous numerical investigations 
appear to indicate that capture of resonance is not strongly dependent upon the 
turbulence model in the cavity [30, 34, 36]. 


2.3 Transformation to Curvilinear Coordinates 

In order to adequately resolve the solid boundary/fhiid interaction, it is common to 
transform the governing equations into curvilinear coordinates which can be body- 
conformal. Specifically, the body is constrained to lie at a constant £, r}, or C level. 
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For a stationary grid, this transformation can be expressed as 

T = E S = Z(x,y,z), il = T]{x,y,z), C = C (x,y,z) 

Application of the chain rule of differentiation yields 

d e d d d 

dx~^ z d^ + 71 x dj 1 + <:i dc 

with similar expressions for the partials with respect to y and 2. The inverse trans- 
formation gives 

d d d d 

aj ~ X( a~x + V( di + 

Again, expressions can be found for the t} and £ partials in a like manner. Represented 
in matrix form: 


_ i 


£x V* Cr 


1 

1 

d_ 

dy 

II 

£y Vy Cy 


d_ 

drj 

d_ 

■ dz . 

s. 

1 

S? > 
Cs 

1 

✓ 



T 

and for the inverse transformation, 


d ‘ 


, 


d_ I 



x i Vi z i 


dx 

d_ 




d_ 

dy 


X Tf J/tj z r) 


dy 

d 




d_ 

dC J 

V 

. X C Vi 2 C . 

✓ 

- dz . 


T -i 


Combining the use of T = (T -1 ) -1 and finite volume metrics, such as those described 
by Vinokur [61], leads to a scheme which is freestream-preserving because of the 
telescoping property. Hence, if the surface normals to a constant f , 77, or £ plane are 
defined respectively as 


s , + i 


+ 5 »,<+jj + 5 «.i+i k 

^(r 7 - r 4 ) x (r 8 - r 3 ) 



CHAPTER 2. NUMERICAL METHOD 


17 


s ;+| = s *j+J i + 5 i/j+^ + s *J+i k 

= ^(r 7 - r 2 ) x (r 3 - r 6 ) 

S k+\ = S z.Jfc+li + S y,it+lj + S z.*+± k 

= ^(r 6 - r 8 ) x (r 5 - r 7 ) 
where the index convention is shown in Fig. 3. 



Figure 3: Hexahedral cell and stencil 
The metrics can then be formed as 

(x = “ !*<*•») = y S r,,+ ! 

= J( x C z r) ~ x v z () ~ y s y,«+£ 
6 = J(x n y< - Xtfr,) = ^,,+1 

Vz = J(yt z t~yt z d = ^ 5 xj+i 


*ly = A x t z <-x<;Zz) = ys vj+ i 

?lz = J( x <y( X (V() 

Cx = J{y^ z r} ~ y^ z t) — 
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Cy J{x t)2( x Z z t]) — yS y j. + i 

G = J^Vr,- X^) = 

These metrics represent the projections of the cell face normal into (x,y,z) space. 
The faces of the hexahedron exactly enclose the discrete control volume, i.e., no gaps 
are permitted at the edges. 

Finally, the Jacobian of the coordinate transformation is equivalent to the inverse 
of the volume, as related by 

_ d{x,y,z) 

J d(CvX) 

= Xt(y n z< — y^z v ) — x^y^z^ — y$Z() + x^y^z^ — yr,Z{) 

= v = g(s,_i +s ; _i + s*_i) • (r 7 - ti) 

Utilizing these metrics in the application of the chain rule to Eq. (1) and subse- 
quent simplification yields 

Q' r + E'AF^ + G' c = 0 

where 

Q’ = QV 

&N ’S = {FnsZx + ENsiy + <-*.VsG) V 
= [EffsS x + FtfsSy + GVs 5 r]t 
Fits = (Fns'Hx + FtfsVy + GnsVz) V 
= + FyvsSj, + C?7vs s i]i 

G'xs — (EnsCx + FtfsCy + GxsCz)V 
= [EtfS s x + Ftfs s y + GtfS s z]k 


Separating the inviscid and viscous portions of the flux vectors, then in the f direction 
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E' ss = E' A E' v , where 

pU 

puU A 

E' = V put/ + £ y P 
pwU A £ : p 
(e + p)U 

Here the contravariant velocity component in £ is U = u^ + u^y + tu^-, without metric 
normalization. The viscous flux can be represented as 


T Tyz^y "b ^~zx^z 
E v V T x y£r + Tyy^y *b T zyt,z 

T~xz£x "b Tyz^y "b Tzz^z 

. (ue' 2 + ue' 3 A we\) + (q x £ x A q y Z y + q 2 C) . 
where the viscous stress terms are evaluated by again invoking the chain rule, and 
the flux in the rj and ( directions are found similarly. The results presented herein 
are implemented using either the thin-layer or the full viscous term treatment. 

The widespread use of the thin-layer approximation, first implemented by 


Steger [62], can be justified from either physical or algorithmic arguments. Physi- 
cally, the neglect of all diffusion processes parallel to the body is similar to that used 
in boundary-layer theory, albeit not as restrictive. Hence, w’hen the viscous effects 
are confined to thin regions along a constant £, r ?, or ( plane, this assumption is valid. 
Regarding the algorithmic argument, the banded matrix structure used in multidi- 
mensional algorithms which sequentially solve a set of unidirectional problems can 
include only these thin-layer terms implicitly. This thin-layer flux in the rj direction, 
assumed to be the body normal coordinate, is expressed as: 

0 

miu v + m 4 v v -f m 5 w v 
m A u n -I- m 6 it^ 

Tn^Urj + m 6 v n -I- m 3 w v 

miuurj A m 2 v v n + m^ww^ -I- m A (uv v A vu n ) + m 5 (uu.^ + wu^) 

-\ -TTleivWr, A WVr,) A K€ p {rj \ A rfy 4- Tf]) 
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where the ( ) denotes an arithmetic mean value and 


3 

II 

■fc 

(& + $ + *) 

, m 4 = 

P 

3** 

m 2 = p | 

[vl + \nl + rjfj 

, m 5 = 

P 

m 3 = p | 

[pi + vl + 

, m 6 = 

P 

3^* 


These viscous flux terms may be found for the remaining spatial coordinates as well. 
The results presented here are implemented using either the thin-layer or the full 
viscous term treatment, as required by the the flow physics. 

Nondimensionalization 

The governing equations may be nondimensionalized by the choice of a length scale, 
denoted by L, and reference values of p, u, and p such as 

Pref Pooi Href \Jpoo I Poo , Pref ~ Poo 

The nondimensionalized variables follow: 

P = PI Pref i P = PI Prof, e = e/{p ref u 2 ref ) 
u = u/u ref , V = V /u re f, W = w/u ref 
t = tu re f/ L, T — p/p, p = Pi pref 

The Reynolds number resulting from this procedure is Re = p re fLu re f /p re f, where 
the (~) denotes a dimensional quantity, and ( ) QO denotes the freestream conditions. 


2.4 Upwind Schemes 


The numerical scheme will be described using first-order terms, following which the 
higher-order extensions will be outlined. The scheme expressed for a cell which has 
a mean flux value on each of the six sides is 


d_ 

dr 





)d(drj 
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+ 

+ 



- f'lj-ij'MdC 
~ k _i)dr)d£ = 0 


In fully discrete form, after dropping the primes for convenience, the governing equa- 
tions can be written as 


+ (F n + \ - F n+1 ^ 


(MM -ML })} = « 


where n denotes the time level in this implicit representation, and A{, Aq, and AC 
are set to unity for convenience. 


These flux terms may be evaluated using a technique which may be broadly classed 
as erther central or upwind. The latter technique is chosen for this study lor the de- 
suable nun.er.cal properties, such as diagonal dominance of the flux Jacobian, and for 
the physical dependence on zones of influence which are inherent in upwind schemes 
Upwind schemes bias the derivative evaluations required to determine the flux 
across fluid cells according to the sign of the characteristic speeds. In this manner 
t ese methods bring the physics of the hyperbolic system, the unsteady Euler equa- 
tions, into the numerical solution process. To facilitate the implementation of these 
upwind schemes, the eigensystem is determined. The similarity transformation which 
diagonalizes the unsteady, inviscid, gas-dynamic equations, shown by Warming et 
al. [63], is outlined as follows 8 ' 


where the rows of T~ l 


dE 

dQ 


= TAT - 1 


are the eigenvectors and 


A - A + + A = drag [£, U, U, U + c, U - c] ||s|| 

using normalized contravariant velocity components. The eigenvalues can be split 
according to, among other splittings, their signs: 

A± = A ± |A| 

2 
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where A is an element of A. 

Two upwind schemes are implemented here to compare the results which may 
be obtained wdth either of the techniques. The initial portion of this discussion will 
be presented unidimensionally for simplicity; the multidimensional extension will be 
outlined towards the end of this section. 


2.4.1 Flux Vector Splitting 

The shock-capturing scheme developed by Steger and Warming [64] revisited the 
classical characteristic procedures. They found that the Euler equations possessed the 
property of homogeneity of degree one for the equation of state used here, meaning 
E(aQ ) = aE(Q). For a vector with this property E = AQ, where A is the flux 
Jacobian given by dE/dQ. Consequently, the flux vector can be split into two parts, 
each physically corresponding to the right and left moving waves. This technique 
resulted in the flux being represented as a combination of the subspaces associated 
with the positive and negative eigenvalues, expressed as 

E = T(A + + A~)T~ X Q = (A + + A-)Q 
= E + + E~ 


where T and T 1 are the right and left eigenvectors of the flux Jacobian matrix .4, 
respectively. The flux across a cell face can be determined by 





“ A t+$Q' + A i+ 


Because the Jacobian at i + 5 is dependent on two states, this solution method now 
diverges from the original Steger- Warming flux vector splitting. The treatment of 
this Jacobian is shown in a following section. Linearization in time can be performed 
in one dimension as follows, extension to the multidimensions is straightforward and 
is omitted for brevity. At the n + 1 time level, 


E ::i = 


or 1 + ^-)ziq?x 

- + £," +i 
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where the implicit change in the dependent variables is given by 6Q n = Q n + l — Q n . 

Note that the Jacobian matrices are frozen at time level n. The remaining flux. E"^, 1 . 

* ^ 

may be obtained similarly. 

To assess the effect on stability of this type of linearization, a procedure developed 
by Barth [65] is applied to this method. Using the semidiscrete form dQ,/dt = — /?,, 
then 


R, = ±(e, h -e h ) 


Ax 

Using frozen Jacobian matrices, the method can be linearized as follows: 



(4) 


where for the first-order subset the Jacobian is a block tridiagonal matrix. The blocks 
along the i th row are 


OR, 

dQi - 1 

dRi 
dQi 
dlL 
dQ ,+ 1 



This scheme is inherently conservative in space because of the telescoping property of 
the finite- volume formulation; analysis of this scheme reveals that it is also conserva- 
tive in time, possibly allowing the use of large Courant numbers [65]. A demonstration 
of this analysis proceeds by writing the scheme as 



6Q = -R n = —A n Q n 


hence 

Q n+1 + AtA n Q n+1 = QT (5) 


In order for the scheme to be conservative in time over a periodic domain, the global 
average of a solution must remain constant for all time, i.e., J2j=i Qt = 52i=i Q " = 
TE[=iQ^ +X - Hence, when Eq. (5) is summed across the domain, the result is that 
the columns of A n must sum to zero. For example, summation across a three-point 
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domain yields 


\Qx + Q 2 + Q 3 ] n+1 a At [/, /, /] 


' B\ Cl . 

n 

' Qi ' 

n+1 

a; b; c; 


Q 2 

= [Q\ A Q 2 A <?3]" 

. ■ ^3 b; 


.Qs. 



where the elements of the middle column of A n are C," = A ~ to , Bl = At,, - A~ 

n 1 o/z 1 l 5/2 3/ 

A 3 = ~^ 5 / 2 i which sum to zero. 


2.4.2 Flux Difference Splitting 

Flux difference splitting methods are based on the Riemann problem, solved exactly 
by Godunov in 1959 [66]. The Riemann problem is composed of m + 1 piecewise 
constant states separated by m wave families. The waves include shocks, contact 
surfaces, and rarefaction fans. For each of the Riemann problem cells, the transition 
of the dependent variables is a function of a parameter family. The solution can be 
found once these transition states are known. Approximate Riemann solvers simplify 
the numerics of the problem by eliminating the iterative process required to find the 
intermediate states. 


Expansion wave 
u-c 



Contact surface 
u 

Compression wave 
u+c 


Figure 4. Riemann problem schematic 


Figure 4 shows a schematic of the Riemann problem with the piecewise constant 
states separated by the appropriate wave families. The flux through the cell face is 



1 

2Ai 

1 

2Ax 


Ei + E i+ i 
( Ei A Ei + 1 


- 

-A|£|, +i ) 


) 
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where |£| = |.4|Q = (.4 + — .4 )Q = [7~(A + — A )T 'jQ. The flux differences associ- 
ated with the + and — traveling waves are 

A£, + + , = (TA+r- 1 ),+!«?.+, -Q.) 

A E-j = (TA-T-') j+i (Q i+ .-C,) 

2 2 

Again utilizing the semidiscrete form dQi/dt = — then 

Ri = ^ [(£<+. - Ei-l) - (A|£| j+ J + A|£|,..J)] 

This method can be linearized using a procedure similar to that described previously 
in Eq. (4). The blocks of the i th row are expressed as 

dR, 1 / dA|£|,_i\ 

dQi-i ~ 2Ax{ - 1+ dQi-x ) 


dRi if dA\E\ t+ i dA\h\, 

dQi ~ 2Ax V dQi ' + dQ, 

dR, __ If 5A|£| i+ i\ 
dQ, +l 2Ax V ,+1 dQ, +l ) 
Substitution of the flux difference splitting expression yields 
dA|£|, + i 0 r. 


dA\E\ 


dQ ,+ 1 


|A| I+ i(<2i+i — Qi 


8A\A\ l+ i 

2 c'v.+i 

These true Jacobians are expensive to compute, and the simplification to approximate 


Jacobians is made as 


dQi+i 


Ml.+i 


Utilizing these approximate Jacobians, the linearization proceeds as 

E 'A = 5^(«: +1 +^i+i 1 -ai£i:4) 

= N. [ ( - 4: + + (- 4 " + > - i- 4 i" + i >«M + 

where E n+1 = E n - f A n SQ and Q n+1 = Q n + (dQ/dt) n At. The fluxes through the 
remaining faces are determined similarly. This scheme can also be shown to obey the 
criterion for conservation in time. 
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2.4.3 Roe Averaging 

In order to determine the Jacobian at the cell face i + some function, A = 
A (Ql , Qr ), must be assumed, where the subscripts indicate left and right states. The 
location of this flux evaluation is one of the differences between finite-difference and 
finite- volume schemes. The evaluation used here is attributable to Roe [67], which 
provides an approximate solution to the Riemann problem. This Jacobian is cre- 
ated through the use of a parameter vector composed of a geometric-like mean of the 
states. The more obvious arithmetic Jacobian forms, such as A = ^(Ai + A R ) or A = 
A-i-iiQL + Qr)) 1 are not conservative forms. Conservative Jacobian forms satisfy 
a (Ql, Qr){Ql - Qr) — El - E r . Stated explicitly, the Roe averaging operation is 

P — \[PlPr 

- _ UiLy/pL + U iR Jp^ 

tt| — 

y/pi + \fPR 

_ Pl u \l + p{Uii + Ujn) + PRUjJl 
Pl + 2p + Pr 

_ Piihr )l + PjSJ}r)L + (hj^jt) + PRjhj^R 
Pi + 2p + p R 

where a (") denotes a Roe averaged quantity, and the latter forms are presented as 
inexpensive alternative expressions. Substitution of (?) for L and (?' + 1) for R allows 
the evaluation of the Jacobian at the intermediary cell face. For the flux vector 
splitting case described earlier, MacCormack [5/] has found this average helps to 
alleviate excessive numerical dissipation in regions dominated by viscous effects. Roe 
averaged values are utilized throughout the development presented here. 

2.4.4 Higher-Order Extensions 

Spatially first-order methods frequently provide inadequate resolution of the flowdield. 
However, the methods discussed above can be extended to higher-order spatial accu- 
racy by modification of the right hand side. In order to assist in the preservation of 
well-behaved solutions near the discontinuities admitted by the strong conservation 
law form of the Euler equations, a total variation diminishing technique is imple- 
mented. If the total variation of a solution is defined as 7V(ti) = '£°1_ 00 |u l+1 - u,|, 
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then a solution which follows 7’V'(u" +1 ) < TV(u n ) is TYD. The TYD constraint can 
be shown to result in diagonal dominance, allowing the use of relaxation schemes. 
In this manner the scheme may be extended to higher space accuracy throughout 
the smoothly varying regions of the field, reducing the accuracy in localities of high- 
gradient and extrema in order to obtain sharp and oscillation-free resolution. These 
methods are rigorously applicable only to scalar nonlinear equations or a system of 
linear equations in one spatial dimension. Application of these schemes to multi- 
dimensional systems of nonlinear equations are generally not TVD. Moreover, it is 
not clear that the higher-order accuracy of the unidimensional problem is retained in 
multidimensional cases. However, the results which can be obtained demonstrate the 
usefulness of the technique. 

Of the several methods which fall into the TVD domain [9], the technique im- 
plemented here is one attributable to Chakravarthy and Osher [68], the development 
of which follows for completeness. In this formulation, the higher-order flux can be 
expressed as a sum of a first-order flux, denoted E t+ 1 , and a flux correction term. The 
flux correction terms are determined by first computing the flux differences across the 
m wave families mentioned previously. Subsequent limiting of these flux differences 
and summation across the wave families results in the higher-order flux. This flux is 
expressed as 


E j+ i = E i+i — — 

3+2 '+2 4 

m ‘ 

. t 

(1 + 40 

4 

™ I 

» S 

(1 + 40 

4 

m . 

. t 

+ (i-40 

4 

' m + 

Z dE H 

. i 


where (~) and (~) indicate a quantity that has been limited,,; is the index denoting 
the wave family, and i is the index assigned to a cell center. Using the notation of P 
for the rows of the left eigenvector matrix, T~ l , and r J for the columns of the right 
eigenvector matrix, T , then the measure of the change in the dependent variables is 

o} + j = ■ (Qi + 1 - Qi) 

The measure of the change in the flux is defined as 


a } = A V = (A j+ + \ ] ~)oP 
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the eigenvalues being split as shown previously. The limited counterparts of these 
values are obtained as: 

a~ 3 = minmod tr~ lt 0o~ t 

o, + i = minmod 1 

= minmod 

■i-j [ «+j i-jJ 

ffili = minmod <r + x , 3 a + , 

This limiter returns the argument of smaller magnitude when the signs are equal, and 
returns zero when the arguments are of opposite sign. This procedure effectively adds 
dissipation locally in regions of high flux gradient and at inflection points. In this 
manner, monotonicity is preserved by preventing the creation of new extrema while 
preserving the global accuracy of the solution. While formal accuracy estimates are 
difficult to ascertain because of the nonlinear application of limiting to different wave 
families, numerical experiments have demonstrated that the global accuracy of the 
underlying scheme is preserved [69]. 

The compression parameter, /?, is restricted according to 1 < /? < and t h e 
limiting operator is given as 

minmod(x, y) = sign(x)(max{0,min[|x|,t/sign(x)]}) 

The compression parameter reduces the amount of dissipation added, the range being 
bounded by accuracy and TVD constraints. Finally, the limited flux difference values 
are expressed as 



“ ^+§ r «+§ 



dE i+\ 

= 

dE’.A 

H 

M|~ 

+ 


=j + 

= ^,_ir _i 

* 2 * 2 

dE i i 

1 2 


This asymmetric limiter is designed to modify the fluxes only in the rapidly varying 
portions of the flow, where nonphysical oscillations are likely to occur. Since these 
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high-gradient regions are confined to thin regions, the dominant solution domain is 
differenced in accordance with the underlying scheme. Variances in the value of the 
compression parameter allow the fluxes to be limited for different gradient levels. This 
implies that use of P max will cause the limiting action to be taken only in the high- 
gradient regions, and lower values of p will result in limiting for commensurately lower 
flux gradients. The variety of schemes which can be obtained using this technique 
are shown in Table 1. 


0 

Unlimited Scheme 

Pmax 

2 nd order TE 

-1 

1 

3 

0 

1 

3 

1 

2 

1 

Fully upwind 

Nameless 

Fromm's 

3 rd Order 

Low TE 2 nd Order 

Central 

2 

5 

2 

3 

4 

5 

oc 

j(A xff xxz 

g(Ax) 2 /nx 

A(Ax) 2 /„ r 

0 

-A(Ax) 2 / IIX 
-I(A xff xxx 


Table 1: Summary of schemes 


Here TE = (J - 0)(Ax) 2 /rrr/4 defines the leading term of the truncation error 
for the unlimited form of the schemes. Local metrics have been used in the above 
method to maintain reasonable computational efficiency, a satisfactory' approximation 
for grids which do not contain rapid variations. 

2.4.5 Viscous Terms 

The viscous terms are treated through central differencing about the cell faces. The 
explicit terms are conventionally differenced after chain-rule expansion, inclusive of 
the cross terms if these diffusion processes are significant for the problem at hand. The 
left hand side does not include these cross terms, and the resultant viscous Jacobian, 
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employing V = [p, u, v, w, e] r as the primitive variable vector, is 




0 

0 

0 

0 

0 

X, dV 


0 

777 2 2 

l~t S £ S y 1 3 

ps x s z /3 

0 

M= ~dt 

TL 

0 

m 2 3 

m 33 

PSyS Z / 3 

0 


0 

777 2 4 

77734 

77744 

0 



0 

"752 

"753 

77754 

™55 . 


m 2 2 = p(^s 2 + + s 2 ) , m 33 = p{s 2 z + |*J + «*) 

4 

m 44 = p(s 2 + s* + -s 2 ) , m 55 = «( s 2 + s 2 y + s 2 ) 


m 52 = U 77122 + W77723 + 17777724 

m 53 = U 777 32 + 7777733 + 7/7 777 34 

77754 = 7777742 + 7777743 + 77777744 


The viscous flux through a cell wall at j + \ is of the form M j+ i(N j+ iQ j+l - NjQj) 
where 



1 

P 


P 0 0 0 0 

-u 10 0 0 

-77 0 10 0 

—w 0010 

— e + |(77 2 + 77 2 + 7i7 2 ) —77 —77 —777 1 


Now, using three-dimensional indices ( i,j,k ), expansion of the block structure gives 
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+ 

+ 

+ 


+ B tj+C k 


At 

V 

At 

~v 

At 


+ ci-..., -C~ , , 

C i ^l c+ Lj^QiJ,k+l 


'V A i+U,k \ 6 Qi+lJ* = ^Qij.k 


where only the thin-layer terms in 77 are shown here. 


2.4.6 Factorization 


The extension of the techniques given above is accomplished through dimensional 
splitting. The method used here is that of Yanenko [70], where the factors are chosen 
in the £, 77, and £ directions. Expressing the three-dimensional equations in compact 
notation as 


At At At 

( I+ - ( A+ r v B+ rC )6Q = AQ 


then the factorization procedure yields 


(/ + ^)(/ + |Ib)(/ + ^cw = aq 


This system can be solved sequentially through the use of intermediary steps without 
loss of time accuracy. Although alternating direction implicit schemes of this type 
offer advantages of vectorization, the system is solved as a sequence of unidimensional 
problems, hence limiting the size of the time step due to stability restrictions [71]. 
The use of this technique here is justified by the requirement of adequate time history 
flow resolution, thus imposing an additional constraint on the maximum time step. 
Application of a line Gauss-Seidel method to the starting shock tunnel test case, 
discussed in Section 5.1.2, confirmed this hypothesis. This relaxation method offered 
a slightly increased stability range, but not enough to offset the additional expense 
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caused by short vector lengths for the shock tunnel problem. Additionally, since for 
the factored scheme the flux exchange occurs at the same time level, the technique is 
conservative, even when convergence at the subiteration level is not attained for each 
time step. 

The expense of both the upwind algorithms is relatively high: 86/is per cell per 
iteration using a single processor on the Ames Research Center CCF Cray Y-MP/832. 
These vectorized codes have computation rates of approximately 140 MFLOPS. In 
addition, the memory requirement is 40 words per cell. Decreased processor times 
may be achieved by many methods. For example, freezing the flux Jacobians for 
several subiterations will offer a processing time reduction of 15% per subiteration, 
albeit at the expense of memory. It is clear that the expense of these upwind methods 
is warranted only for problem classes in which the improved resolution is critical. 


2.4.7 Newton Iterative Technique 


Reduction of the linearization and factorization errors is achieved by a Newton iter- 
ative method of the type described by Rai and Chakravarthy [74] and Rogers and 
Kwak [75], albeit with the addition of allowance for a varying step size [14]. Assuming 
that the initial guess lies within the radius of convergence, the right hand side is con- 
verged to an arbitrary' accuracy while holding time fixed. Since the right hand side 
includes the higher-order difference representations of the Navier-Stokes equations, 
linearization and factorization errors are eliminated at convergence. The method is 
discussed below where m is the Newton iteration index and n is the conventional 
index denoting time level. Discretizing Q t + E z = 0 gives 


1 (Qn-U.m+l _ gn + l,m^ _ ^gn+l,m+l _ gn,m + 1 _ ^n+l,m _ gn^j 


gn+l,m+l l_(gn+l,m _ gnj 

At 

where the solution is converged at time level n, hence <2 n,m+1 = Q n . Defining 6Q' = 

gn + l,m+l _ gn + l.m^ then 


1 


I6Q' 


= atq; 


n + l,m+l 


- (Q 


n+l,m 


+ l _ ^Qn + \,m 


Q-) 

-<?") 
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Linearization at iteration level m + 1 gives 


^n + ^m + l ^n + l*m 


+ 


'dE z 

dQ 


n+ 1 ,m 


6Q' 


where the flux Jacobian has been frozen at iteration m. Substitution yields 

rs 

I6Q' = — Ar£ , ” +1,m - AT-^-A n+hm 6Q' - (Q n+l ' m - Q n ) 

ox 


Rearranging results in 

bQ' = -(Q n+1 ’ m -Q n + Ar£’” +1,m ) 

which reverts to the standard noniterative form + ^.A n j bQ 1 = — when no 
subiterations are taken, as can be seen by substitution of n for n + 1, m. 

The temporally second-order accurate representation is found by extension of the 
above procedure. Using a three-point backward time stencil derived from a standard 
Taylor series approach, 


I + 

Ax 


Qt = C 0 Q n+l + C x Q n + CaQ"- 1 


where 


Co 

C 2 


l — o 

(1 - o)At 2 + AT! 
-1 


, c, = 


(1 — o)At 2 -f- A Tj 


(1 - o)At 2 + Ar! 

and o = (1 A Arj/ A t 2 ) 2 . The elapsed time between the n — 1 and n time levels is 
given by Ati and between n and n + 1 is A r 2 . Rewriting at iteration level m + 1, 


C 0 (Q 


n + l,77i + l 


Q n+1 ' m ) = Qj 1+1 * m+1 _ (C 0 Q n+l ’ m + C,Q n + C 2 Q n ~ l ) 


Finally, 


1 + 


1 


CqAx 


A 


n + l*m 


bQ' = -(Q n+1 ’ m + 0-Q n + - -^£ r n+1 ’ m 

Co Co Co 


which reduces to 


- + bQ ' = 2 aT (Q n — Q n 1 ) — E* for the case of no iterations 

with fixed time step size. The formulation given above allows the use of a time step 


_3 1 
2A 
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size which is a function of time, but is fixed at each step for the entire domain. The use 
of a variable time step size allows the solution to progress using a constant Courant 
number, possibly preventing inadvertent divergences. Higher-order accuracy in time 
may be obtained by extension of the above technique, albeit with additional memory 
requirements. 

The assertion that this technique reduces the factorization and linearization errors 
is substantiated as follows. The right hand side of the method contains the discretized 
governing equations in their pure form, that is, without the numerical approximations 
utilized to attain rapid convergence. The left hand side allows the use of large time 
steps by relieving the Courant-Friedrichs-Lewy stability constraint. Deferring the 
question of uniqueness, if a set of dependent variables is found such that the right 
hand side is satisfied, then this field is a solution to the discretized equations regardless 
of the approximations made to arrive at that set. 


2.5 Diagonal Scheme 

Resolution of the transient wave-field of the blast- w'ave/target interaction problem 
class benefits from the use of upwind methods. In contrast, for the SOFIA effort a 
combination of the low transonic regime, the complex geometry, and the large prob- 
lem size required an efficient integration scheme. The algorithms used for the SOFIA 
effort, coded by Buning and Chan [58], are implemented w'ithin the Chimera over- 
set grid framework [72]. The solutions were obtained using a diagonal scheme [73], 
using spatially varying time steps for steady state computations, and fixed step size 
for unsteady flow simulations. The code utilizes the conventional dependent variable 
vector, Q = [p y pu,pv, pw,e] T , and pseudo-finite-volume metrics. Euler implicit time 
marching and second-order central spatial differencing were used for the computa- 
tions presented here. Computations were performed on the Numerical Aerodynamic 
Simulator (NAS) Cray Y-MP/832 using SSD, at an expense of lAps per point per 
iteration. 
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2.6 Geometry Treatment 

Geometric modelling and grid generation is a significant portion of the effort spent 
in obtaining the flowfield about any reasonably complex geometry. A structured 
approach is utilized for this study, with the body-conformal internal grids generated 
using the elliptic techniques of Thompson, et al. [76], Thomas and Middlecoff [77], and 
Steger and Sorenson [78]. The external flow domain was discretized via algebraic and 
hyperbolic means, the topology w r as chosen to allow the use of these grid generation 
methods. The grids used in this investigation were generated using codes written 
by Steinbrenner, Chawner, and Fouts [79], Chan and Steger [80], and Atw r ood and 
Vogel [81]. A discussion of the treatment of the surface, the grid topology, and the 
grid strategy is given below. 

Surface Modelling 

The geometry used for the SOFIA configurations utilized clipped wings to emulate 
the geometry used in a specially designed experiment [16]. The use of clipped wings 
in the wind tunnel test allowed a cavity of more realistic size to be studied. The 
fuselage, wing, fairing, nacelles, and telescope geometry were obtained from CAD 
databases. Positioning errors in the database were corrected using blueprints. 

The process of generating the more complicated grids, e.g., the quiet SOFIA 
configuration with telescope (configuration 100), warrants additional comment. Ex- 
tensive wind tunnel testing [16] resulted in a hand-formed ramp and aperture, shown 
in Fig. 5a. This geometry was subsequently laser digitized [82], resulting in a data set 
of the form shown in Fig. 5b. These data, accurate to approximately 0.2 mm (0.6% 
of cavity length), were then converted into a form suitable for the surface grid using 
a standard CAD package [83]. 

Surface definition via bicubic surfaces in regions of high curvature can cause local 
oscillatory behavior [79, 81]. Along overlapping surfaces this property is manifested as 
C°, or jump discontinuities at zone boundaries. The problem is ameliorated through 
bilinear projection from one zone boundary to another [84]. The distance of projection 
is typically five orders of magnitude less than a characteristic geometry length. 
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(c) 

Figure 5: Geometry acquisition: (a) model, (b) laser digitized configuration, and (c) grids 
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Topology 

The overset grid topology scheme was chosen for its geometric flexibility and its ability 
to allow refinement of individual zones. The use of a different grid for each component 
of the geometry simplifies changes that will occur as the design matures. In fact, the 
geometries with cavities were built upon the clean configuration grids, providing a 
savings of many man-hours. The topology' w r as chosen to allow rapid evaluation of 
new configurations and permit simple specification of turbulent wall and shear layer 
regions. 

Grids 

The SOFIA configuration y + values for the first grid point away from the wall were 
generally about 4.0, the farfield boundary was placed at 20 fuselage diameters, and 
the outflow 10 diameters downstream. Damping of acoustic waves at the farfield 
boundary was achieved by the use of large cells which were unable to support the 
high frequency waves. 

The clean SOFIA 747 configuration, without a cavity, was modelled using four 
grids for the half-body: one each for the fuselage, wing, wing tip, and nacelle. The grid 
point count was approximately 4 x 10 5 . The fuselage grid was refined in anticipation 
of the cavity to provide similarly sized cells in interpolation regions. 

The untreated aperture geometry, configuration 25 of the wind-tunnel test, was 
gridded by reflecting the four grid zones described above and adding two for the 
cavity. The term untreated refers to the lack of geometry modifications which can 
eliminate cavity resonance. The fuselage zonal boundaries were shifted meridionally 
to move interpolation away from the cavity region. The two additional grid zones 
consisted of an outer cavity grid surrounding the cavity' region and an inner cavity 
grid which included the cavity walls and the shear layer region. The outer zone was 
utilized to isolate the cavity unsteadiness from the global solution. The total grid 
point count for this case was about 1.2 x 10 6 distributed in 10 zones. 

The treated aperture geometry, configuration 100 of the wind-tunnel test, was 
modelled by the addition of seven grid zones to the clean case: one each for isolation, 
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aperture wall, shear layer, cavity wall, telescope tub, secondary mirror, and an inner 
ramp grid. The total grid point count for this case was about 1.8 x 10 6 in 15 zones. 


Chapter 3 

Boundary Conditions 


The flow solver block implicit boundary conditions are implemented in a manner con- 
sistent with the flux split linearization described earlier. The inviscid and viscous im- 
permeable wall conditions are prescribed similarly to those given by MacCormack [ 85 ]. 
Although the following procedures are presented for a cell face which lies along a con- 
stant 77 plane, the procedure may be generalized for application to any cell boundary. 
Finally, the characteristic inflow and outflow boundaries are discussed. 

The inviscid, impermeable wall boundary condition is described for a pair of cells 
between which the surface lies, depicted in Fig. 6. In the following discussion, the cell 



39 




CHAPTER 3. BOUNDARY CONDITIONS 


40 


above the wall will be denoted by subscript 2, the cell below the wall by subscript 1. 
At the centroid of cell 2 the velocity is expressed as 

V 2 = u 2 i + V2J 4- w 2 k 

and at the cell wall the surface normal is 

^ 1 vail — 71 rr(^x® 4” SyJ 4" 

ll s ll 

(s x i 4” Syj 4" Szk)!^// 

where ||s|| is the vector magnitude. Hence, the velocity component normal to the wall 
is 


K>2 — ||Ki2||S‘ u , a // 

= (us x + vs y + ws z )(s x i 4- sj 4- s z k) wal , 

^ ««* 

Since, V = V t + V n , then the tangential velocity component is 

Vt2 = V 2 -Vn2 

= [u - s x (us x + vSy + u;s-)]i 
4-[v - s y (us x + vs y + tus 2 )]j 
4-[u/ — s z (us x + vs y 4- u/s 2 )]k| 2 


The flow tangency condition is satisfied by V a = V l2 and V ni = -V n2 . Total energy 
and density are found from reflection as even functions through the relation R6Q1 = 
ER6Q 2 , where E = dia^[l, 1, — 1, 1, 1] and the rotation matrix is found from an 
eigenvalue problem, the eigenvectors being the rows of 


R = 


1 0 0 

0 4“ <$ 2 ) S x 

0 S x Sy 

0 s z s z 

0 0 0 


0 0 

s x 0 

s z 0 

—(s x 4- Sy) 0 
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A Riemann problem can then be solved at the wall to determine the flux from Ej=i — 
B~Q 2 + B + Q\. This amounts to the wall being represented as a contact discontinuity 
by constraining the contravariant velocity to vanish. Implicitly, this results in 

6Qi = R~ 1 ER\ wa it6Q 2 ( 6 ) 


where 


R~ 1 ER 


1 0 0 
0 1 - 2s* —2s x s u 

0 -2s y s x 1 - 2 s* 
0 -2s z s x —2s z s y 
0 0 0 


The block tridiagonal system may be written as 


0 0 
— 0 
— 2s y s* 0 
1-2 s* 0 



B[ C[ 

A 2 f ? 2 C '2 



' 6 Q 1 ' 


■ A Q, ' 


SQ 2 

II 

aq 2 


Now the change in flux across an arbitrary cell wall boundary is given by AE wa u = 
A+6Q\ + A~SQ 2 , or the sum of the changes in the flux contribution from the positive 
and negative moving waves. Substitution of Eq. ( 6 ) yields 


A E wa it = (A + R~ 1 ER + A~)6Q 2 


and it can be seen that dependence upon 6Q X has been eliminated. Hence, the block 
tridiagonal system may be represented with embedded boundary conditions as 


' B" C\ 


■ 6 Q 2 ' 


‘ A Q 2 ' 

^3 B' 3 C ' 3 


SQ 3 

— 

aq 3 


where B" is the appropriately modified Jacobian. 

The viscous impermeable wall imposes additional constraints on the specification 
of the wall flux. Again utilizing a primitive variable vector V = [p, u, v, w, e] T , then 


dR, 

<5Vi = diag[l,t u ,t v ,t w ,t]SV 2 = -^rybV 2 

u V 2 
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for a wall face at | . In this form the toggles t u ,t v A Wl and t are set at 1 or — 1 for a slip 
or a no-slip condition, or adiabatic or isothermal wall, respectively. This may be seen 
by simply rearranging expressions of the form u wa u = i(ui + u 2 ) or u wa u = uj = u 2 . 
Having already specified the impermeable wall conditions earlier, only the viscous 
terms at the wall are of present concern. Looking at the terms of the form 

A Q 2 = ^(-MN X 6Q X + MN 2 6Q 2 + • • •) 


then substitution of the ■wall relations above leaves a term 


_ At 

Q 2 - yj 

At 

- V2 


dV x dV 2 dV 2 

-M^AJ-6Q 2 + M—^-6Q 2 + 


M 


dV 2 dQ 2 

(l-W) 

V dV 2 ) 


dV 2 

dQ 2 


dQ 2 

SQ 2 + 


which is subsequently embedded into the block structure. The dependent variables 
within cell 1 are specified according to boundary-layer theory, holding the pressure 
gradient zero normal to the w r all. The remaining variables follow from fluid and 
thermodynamic relations. 

The inflow and outflow boundaries are specified according to characteristic theory 
for generality. Linearization of T~ l Q t + AT~ l Q^ = 0 for the forward differenced inflow 
condition yields 

[/ - = — yAf-HQ",., - Q") 

where the modified eigensystem matrices are computed as 


A = diag [o, 0, 0 , 0 , (1 — t in )(U — c)] ||s|| 


f~ l = 


dpr/dQ 

dT T /dQ 

dv/dQ 

dw/dQ 

(1 — fin )^5 + tindu/dQ 


Here t in is zero or unity for subsonic or supersonic inflow. The specified variables are 
chosen such that a unique set of flow quantities are given at the entrance. 
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The outflow condition is specified in a like manner; however, there are now at least 
four characteristics linking the domain with the boundary. Backward differencing 
about i 4- 1 and specification of static pressure at the exit results in the modified 
matrices 

A = diag [0, U, U, U + c, t out (U - c)] ||s|| 

h 

h 
h 

u 

(1 — tout)dp/dQ + touth 
In the above development, the eigensystem is evaluated at the boundary face in 
question, maintaining consistency with the interior treatment. 

The strongly unsteady blast-w'ave problems investigated here revealed that the 
use of block implicit boundary conditions resulted in significantly enhanced conver- 
gence. This beneficial effect is caused by the faster signal propagation arising from 
the incorporation of the boundary conditions within the linear system. However, for 
the cavity flows explicit boundary condition implementations were used for coding 
simplicity. Comparison of the computed results w’ith experiment show’ satisfactory 
resolution of the moderate unsteadiness present in transonic cavity flows. 

Several of the cavity cases used characteristic boundary conditions holding mass 
flow, total enthalpy, and flow angle constant. Subsonic inflow conditions, for example, 
are related to the interior of the flow by the u — c characteristic: 

dp du ( dp du 

-- P c- = -(u- c )^- pc - 

which may be rewritten as 
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Pi + dp 

m 





Implementation of the boundary conditions, unless otherwise noted, are as follows: 
viscous impermeable wall conditions are no-slip, zero normal pressure gradient, and 
adiabatic; information transfer across overset mesh boundaries is implemented using 
non-conservative trilinear interpolation of Q. Treatment of the farfield boundaries is 
case dependent and is noted in the results section. 


Chapter 4 

Geometrical Aero-Optics 


The objective of the numerical simulation of the flow about the SOFIA airborne ob- 
servatory is to design a safe configuration which will have the least detrimental effects 
upon the optics. Towards this goal, the following transonic cavity flow problems were 
divided into three sections. First, the unsteady interaction of the external flow with 
the cavity requires time-dependent solutions to the Revnolds-averaged Navier-Stokes 
equations. Second, the shear layer growth rate is strongly dependent on turbulence 
effects which must be modelled due to the grid coarseness. The final portion of the 
problem is the application of the optical model to the unsteady density field in the 
shear layer to determine seeing quality. 

The variation of the speed of light through gases is primarily a function of the 
density field. This fact has been extensively used to benefit the study of fluid physics, 
as exemplified by use of schlieren, shadowgraph, and interferometry techniques. How- 
ever, the objective of the present effort is to quantify the wavefront distortion of a 
beam of light propagating through the shear layer. This distortion is computed using 
the history of the density variations within the shear layer to predict fluctuations in 
the optical path length via geometric optics. This will in turn allow prediction of the 
telescope resolution limits due to seeing and thus contribute to the telescope design 
specifications. 

The geometric optics model developed here assumes that the impact of the fluid 
density on the optical field may be computed by casting light rays through a field 
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Figure 7: Geometric optics: (a) partitioning of hexahedrons into tetrahedrons, and (b) 
intersection and refraction procedure 

discretized into tetrahedrons. Diffraction effects, which become important when the 
wavelength approaches the scale size, are neglected. The simplifications afforded by 
the use of planar facets and piecewise continuous media are utilized by tesselating, 
or partitioning, each hexahedron of the flowfield into five tetrahedra as shown in 
Fig. 7. It should be noted that other tesselations were found to be more robust for 
thin warped cells [86], however for the shear layer grids used here the five tetrahedra 
decomposition is well-behaved. 

Application of the geometric optics code to two preliminary test cases was under- 
taken to determine sensitivity of the optics code to the above non-unique tesselation. 
The parallel emergence of the rays after propagation through a plate and a prism of 
index of refraction n = 2.4 suggests that the results are relatively insensitive to the 
method of tesselation. 

The problem can now be divided into three steps: 1) propagation of the ray, 2) 
intersection of the ray with a facet, and 3) refraction. Solution for the point contained 
in both the facet, p, and the ray, q,, expressed parametrically as 

q i(t) = d + et ; p(u, w) = a + bu + cxv 

results in the intersection in parametric coordinates, which is found from the dot 
product of the normal with the ray and the surface: 

(b x c) • p = (b x c) • q { 

_ (b x c) • a — (b x c) • d 
(b x c) • e 
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(c x e) • d — (c x e) • a 

u = 

(c x e) • b 

(b x e) • d — (b x e) ■ a 
U (b x e) • c 

where the vector coefficients are found from boundary conditions: 

a = p(0, 0) 
b = p(l,0) — p(0,0) 
c = p(0, 1) — p(0,0) 
d = q.lO) 
e = q,(l)-qi(0) 

Specification of the light ray origin and a direction initializes the problem. Following 
the search for the initial hexahedral cell in which the ray originates, the tetrahedron 
within this hexahedron must be computed. First, the shortest intersection distance 
of the ray with the 16 planes which compose the hexahedron is computed. Then 
the dot product of the ray with the fourth vertex of the closest plane determines the 
origin tetrahedron. Subsequent intersection and refraction processes are a marching 
procedure. The optical path length ( OPL ) is found from 

OPL = J n(s)ds « n, As, 

j 

where n(s) is the index of refraction as a function of position along the ray, s. The 
variation of the OPL over the aperture gives a measure of the wavefront error caused 
by the shear layer. 

The refraction process is determined according to Snell’s law as shown in Fig. 7, 
where the planar interface, p (u,w), separating the media and the incident light ray, 
q,-(f) [87, 81] are depicted. Generalization to three-dimensions is accomplished by 
rotation to the osculating plane, which includes the surface normal and both the 
incident and refracted rays. In this osculating plane, a rotated local coordinate system 
is defined: 

q. = |q.„l fi + |q.<|t 
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where 

lq.nl = |q, | COS 6>. and fq. J = |q,| sin <9, 

- _ q, — n cos 9 t 

sin Qi 

Application of Snell’s law nisin#, = n 2 sin0 r where n = 1 + Q- 2 — results in an 

PSTP 

expression for the refracted ray: 

q r = n cos 0 r + t sin 6 r 

The local index of refraction, rij, is found by arithmetically averaging the densities 
at the four vertices of the tetrahedron, where only one vertex changes as the ray 
propagates to a neighbor tetrahedron. The Gladstone- Dale constant, /?, is a function 
of the media and of the wavelength. Using air as the media and a wavelength of 
A = A d = 5893A, then /? = 2.92 x 10 -4 . The Cauchy formula can also be used: 

0 = [2875.66 + 13.412/( A 2 x 10 -8 ) + 0.3777/(A 4 x 1(T 16 )] x 1CT 7 

The values of 0 used in the present computations were chosen to match those used in 
the reduction of the experimental data. The wavelengths of interest for SOFIA range 
from the near infrared, l^m, to the microwave, 1 mm, where optical distortion can 
be seen to be more severe for shorter wavelengths. 

Finally, to obtain a measure of the loss in irradiance due to the fluctuating density 
field, the OPL for vacuum conditions is subtracted from the OPL through the gas 
to yield the optical path difference ( OPD ). The value of < OPD' > is computed 
using a sequence of OPD's at a fixed station. Using the root-mean-square wavefront 
distortion < OPD' >, the phase distortion (4*) is found from $ = The Strehl 
ratio, given by 



is a measure of the peak intensity to which a beam can be focused. The computational 
expense of the procedure outlined above is currently 250/xs/hexahedron/ray on a 
single processor of the Numerical Aerodynamic Simulator Cray 2. 

In these studies, the effect of fluids upon the optical field is determined through 
prismatic modelling of the density field. Integration of the equations of motion 
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through a trilinearly varying density field could also be implemented. The use of a 
six-tetrahedron decomposition would eliminate the present requirement of a checker- 
board cell arrangement to prevent gaps between hexahedrons. The modelling of the 
wave-like nature of light could also be implemented via the parabolized Helmholtz 
equation if diffraction effects were deemed significant [88]. 



Chapter 5 

Results and Discussion 


5.1 Blast- Wave Results 

The methods introduced in the previous sections are applied to test cases which 
demonstrate the capabilities of the algorithm. The viscous term treatment in a low 
Mach number regime is shown by the Couette flow problems, which are compared 
to Similarity solutions and previously obtained numerical results. Demonstration of 
the inviscid term treatment is shown by capturing of transient discontinuities for a 
shock tunnel start-up problem. The three-dimensional results are compared with an 
experimental study of a hemicylinder mounted in a shock tube. 

5.1.1 Couette Flow 

The Couette flow problem is used to compare the present methods against the method 
of Beam and Warming [7] and the similarity solution as given by Schlichting [89]. The 
results for the two upwind methods fall virtually on top of each other, and n indicates 
the time step. The solutions shown in Fig. 8 were obtained using quiescent initial 
conditions and viscous boundary conditions with no-slip adiabatic walls. Both of 
these cases were implemented in the thin-layer form at a Reynolds number of 6.4, 
based on the distance between the plates, equal to 10~ 5 feet. During the course of 
these solutions, slightly more than an order of magnitude drop in \6p\ max per two 
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subiterations was observed in the (3 x 10) cell domain. The Courant number used 



Figure 8: Couette flow case: (a) impulsively started and (b) oscillating plate 


for the oscillating plate calculation was approximately 10, indicating the viability of 
these types of unsteady computations at Courant numbers greater than unity. The 
Courant number was computed using CFL = ^max[\(U, V, t? r )| + c]||s|| over each 
cell in the domain. Identical results were obtained using both the two- and three- 
dimensional implementations in all directional permutations. Results reveal slightly 
steeper gradients than that of the conventionally differenced scheme or the analytic 
solution, a possible consequence of the handling of the boundary' conditions or the 
viscous term treatment. In addition, this case was found to be insensitive to the 
choice of the higher-order flux correction terms, possibly because of the dominance 
of diffusive effects. 


5.1.2 Shock Tunnel Start-up Problem 

The third test case evaluated the inviscid term treatment through the simulation of 
the transient starting process of a planar shock tunnel. The (300 x 60) cell domain 
is shown in Fig. 9. The solution of the Euler equations is presented in Fig. 10 as 
a comparison of experimental and numerical shadowgraph images, the former due 
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Figure 9: Shock tunnel case: 300 x 60 cell grid 


to Amann [90]. The computed shadowgraph function is proportional to V 2 p and 
thus acts as an amplifier of the density gradient across contacts and shocks, for in- 
stance. This solution was obtained using Roe flux difference splitting with <p — 1/3, 
the upwind biased flux evaluation. The Steger- Warming flux evaluation with Roe 
averaging was found to be moderately less stable, but no significant differences in the 
results were found for this case. The maximum compression parameter was used, and 
the entropy-fix parameter used in Harten’s formulation [9] was set to 0.15. Discon- 
certingly nonphysical solutions were produced for smaller entropy-fix levels, possibly 
associated with an entropy-violating condition. The problem was initialized with a 
moving shock propagating to the right at a Mach number of 2.97, while the boundary 
conditions were specified as impermeable inviscid along the walls and fixed for the 
inlet and exit. For the maximum Courant number of four used here, four subiter- 
ations were chosen per time step based on a subjective judgment of discontinuity 
sharpness. The |<5p| moi was observed to drop approximately an order of magnitude 
over the course of these four subiterations. 

Physically, this nozzle starting process generates a high enthalpy reservoir of more 
than 50 times the initial pressure, while the density increases by 11 times the initial 
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Experiment, Amann (Ref. 90) Euler 


(a) Primary shock reflected into reservoir 

(b) Swallowed primary shock 

(c) Rearward facing shock being swept downstream 

(d) Reflected shock system 

(e) Mach line generated from C 2 discontinuity 


Figure 10: Shock tunnel case: shadowgraph comparison for a Mach 3 planar shock 
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state. This reservoir provides the energy necessary to generate high Mach flows 
dow r nstream of the diverging nozzle region for short durations. The ensuing reflections 
of the shock with the nozzle wall reveals the complexities of the shock-shock and 
shock-contact interaction. In particular, it can be seen that the development of the 
rearward facing shock, which is directed upstream while being swept downstream, 
is resolved. At later times, the finer scale fluid motion between the primary and 
rearward facing shocks is, for the most part, lost because of grid coarseness and 
attendant numerical dissipation. However, increasingly fine structures are captured 
as the grid is refined. 

5.1.3 Shock Tube Blockage Study 

The viability of the technique in three-dimensions is shown by the final test case. 
These results are intended to replicate the conditions in an experimental study of 
a blast-wave encounter with a hemicylinder target in a shock tube by Kingery and 
Bulmash [91]. The experimental test configuration and pressure transducer locations 
are shown in Fig. 11. In order to estimate the costs and benefits of inviscid versus 
viscous simulations, the flow about this geometry was computed using both the Euler 
and Navier-Stokes equations. However, the expense of these three-dimensional simu- 
lations permitted the use of only one of the inviscid flux evaluation methods; the Roe 
flux difference splitting was chosen. 

The simulation was initialized as a translating planar shock before diffraction over 
the cylinder began. Initial conditions were specified as: 

M, = 1.518, Re/m = = 23.3 x 10 6 /m 

T x = 288.17 K, poo = 101.3 x 10 3 N/m 2 

Boundary conditions are specified for the viscous, single zone computation as shown in 
Fig. 12. In the shock-tube direction, the ^-direction, extrapolation is used. This non- 
physical extrapolation is adequate for the duration of the early interaction. However, 
solutions at larger times are suspect, where times after the shocks have propagated 
through the boundaries are defined as large. Additionally, the use of an advancing 
front boundary is enabled because of a priori knowledge of the grid structure and the 
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Figure 12: Hemicylinder case: central region of the 78 x 50 x 25 cell (a) inviscid and (b) 
viscous grids 


primary shock speed. This simple time-dependent boundary reduces the computation 
time by using the fact that nothing occurs ahead of the blast- wave. In the 77 -direction, 
the lower boundary defines the surface geometry of the hemicylinder, and hence is 
specified as a no-slip isothermal wall. The top of the domain in the 77 -direction, corre- 
sponding to the inner radius of the shock tube, is specified as an inviscid wall, based 
on the assumption that the viscous effects on this surface have negligible influence on 
the results. Finally, the C-direction boundaries are treated using the viscous condition 
along the floor of the tube. Symmetry conditions are used along the plane running 
along the longitudinal axis of the cylinder and normal to the floor. To simulate ex- 
perimental conditions, the wall temperature was set equal to the temperature of the 
quiescent flow prior to primary shock arrival. The viscous grid has normal spacing 
of approximately 10 -5 meters at the viscous walls. The Euler compu*->tion used the 
inviscid boundary conditions previously discussed where appropriate. For this Euler 
grid, since the areas of the faces corresponding to the geometric axis singularity are 
zero then F • s is also zero. Compensation for the round-off error inherent in the grid 
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was implemented by eliminating those face areas which fell below a specified toler- 
ance. The grids and boundary conditions for these cases are partly shown in Figs. 12 
and 13. 



Figure 13: Hemicylinder case: symmetry plane of the central part of the viscous grid 

The inviscid computation used p = 1/3, (3 = 4, Harten’s entropy fix parameter of 
10 -4 , Roe flux difference splitting, one subiteration per second-order accurate time 
step, and a Courant number of 15. The solution was obtained in 1500 time steps 
without any change of parameters. 

The viscous computation used the same flux evaluation as above with the addition 
of the second-order accurate full viscous terms. Because of the viscous spacing, the 
Courant number utilized was 1 0 4 , allowing the solution to be obtained in 6800 time 
steps with no subiterations. In contrast to the inviscid simulation, the advancing 
front boundary condition was utilized in this case. 

Results, given in Figs. 14 through 17, show that the primary shock is captured 
over two to three cells, the large physical thickness obtained is an artifact of the 
coarse grid used. Adaptive gridding methods would help maintain sharp shocks, but 
the anticipated expense of these methods precluded their use here. Figures 14 and 15 
show comparisons of the numerical and experimental pressure histories. It is seen that 
the peak overpressure is underpredicted by 10%, possibly owing to the coarseness of 
the grid which in turn thickens the shock. These computed surface pressures were 
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extracted from the domain through the use of a Newton search in three-space for 
the cell in which the given ( x,y,z ) probe coordinate fell [92]. Subsequent trilinear 
interpolation over the cell, where the uniform parametric coordinates {u,v,w) are 
determined from the positions of the vertices of the hexahedral cell, allows the pressure 
to be computed. Inherent in this first-order approximation lies the assumption that 
over a discrete cell the variation of pressure is linear in space. 

Figures 16 and 17 show portions of the viscous simulation at selected times. Phys- 
ically, the interaction process begins with the normal impact of the incident shock 
with the front face of the hemicylinder. At this time, peak overpressures of six times, 
and densities of four times that of the quiescent state are generated along this for- 
ward face. As the shock diffracts over the sharp corner of the target, a separation 
bubble forms, which eventually envelops a large portion of the circumferential face of 
the body. This vortical motion is depicted in Fig. 17 by instantaneous streamlines. 
A supersonic pocket is generated as the air negotiates the sharp corner as it rushes 
from the stagnation region left in the wake of the upstream propagating reflected 
shock. The next significant event occurs as the shock diffracts over the rearward face, 
shedding a strong vortex sheet while an expansion wave propagates away in a pattern 
which grows with time. The diffracted shock then impacts the floor of the shock tube, 
reflecting it upwards, while the shock which diffracted over the circumferential face 
reflects inwards from the outer walls of the tube. A simplified sketch of the interac- 
tion process is shown in Fig. 18. The subsequent diffractions and reflections result in 
the interaction of shocks, expansion fans, vortices, and developing boundary layers. 
From experimental evidence, this gross unsteadiness does not dissipate for more than 
15 milliseconds after the interaction event begins. However, the primary shock passes 
from the test section 5 milliseconds after the initial target interaction; therefore, the 
computation is stopped at that time. 

The effects of the viscous terms are seen by comparing the pressure histories in 
Figs. 14 and 15. While the pressures along the upstream face are largely unchanged, 
the circumferential and downstream faces are significantly affected by viscosity. The 
large separation along these faces causes low pressure regions due to this vortical 
motion. This phenomenon is more accurately captured in the viscous simulation, as 
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Pressure [atmospheres] Mach number 



Figure 16: Viscous hemicylinder case: symmetry plane pressure [atm.] and Mach contours 
at (a) t=l.lms, (b) t=1.7ms, and (c) t=6.2ms 
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Figure 17: Viscous hemicylinder case: instantaneous streamlines at (a) t=l.lms, (b) 
t=1.7ms, and (c) t=6.2ms 
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Figure 18 : Hemicylinder case: schematic of shock interaction 


may be seen by inspection of the pressure histories at probe 11. Differences between 
the experimental and the present results may be due to poor capturing of the vortex 
strength owing to grid coarseness. However, the higher-order behavior of the method 
used here attempts to reduce the need for finely spaced meshes. In addition, the 
occurrence of deformation of the shroud wall is thought to be a possible event during 
the experiment, and could adversely affect the comparison between experiment and 
computation [93]. 

A limited cost/benefit study of the Euler versus Navier-Stokes equations was also 
performed for the hemicylinder case. For approximately 5.5 ms of flow history on a 
(78 x 50 x 25) cell grid, the Euler computation consumed 7.6 processor hours, while 
the viscous simulation required 18.2 hours. From these results, the somewhat more 
accurate solution given by the Navier-Stokes simulation may be worthwhile. This 
is particularly true if the flowfield behavior after the direct interaction the primary 
shock with the target is important. 
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5.2 Cavity Flow Results 

Validation of the diagonalized code was accomplished by evaluation of two- and three- 
dimensional cases related to the transonic aero-window problem. Numerical results 
for free shear layer and rectangular two-dimensional cavity flows were compared with 
analytic and experimental data to evaluate the capability of capturing the fundamen- 
tal physics. Analysis of the SOFIA configuration simulations, including evaluation of 
optical distortion is also presented. 

5.2.1 Free Shear Layer 

A series of numerical experiments was performed using a two-dimensional shear layer 
as the test case. Sensitivities of mean and time-varying quantities to changes in time 
step size, fourth-order dissipation levels, and grid refinement were determined. Addi- 
tionally, partial validation of the algebraic turbulent shear layer model was determined 
through comparison with similarity solutions and experimental data. 

The computational domain for this case includes a two inch long splitter plate 
embedded in a channel, with initial conditions specified as a discontinuous step at 
the channel centerline. Shown to scale in Fig. 19, the channel extends 30 inches 
downstream of the splitter plate trailing edge, and five inches above and below the 
plate. Inviscid walls were specified for three inches upstream of the viscous splitter 
plate and for the channel walls. The inflow and outflow conditions were implemented 
using one-dimensional characteristic relations holding mass flow, total enthalpy, and 
flow angle fixed at the inlets, and fixing pressure at the exit plane. The boundary 
layers on the splitter plate and the shear layer were turbulent. Reynolds number 
based on the mean velocity of the streams and the length of the splitter plate was 
6.7 x 10 5 . 

The results for three grid refinement levels are shown in Fig. 19 along with 
Gortler’s similarity solution. The velocity profiles are taken 10 inches downstream of 
the trailing edge of the plate. The solution can be seen to become grid independent 
when the grid becomes finer than approximately 20 points across the layer. The Mach 
number ratio for this case was 0.2/0. 8 and a — 20.7. Eddy viscosity was observed to 
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Figure 19: Shear layer case: (a) velocity profiles of differing grid resolution compared to 
similarity and (b) variation of spread rate with velocity parameter 

grow linearly in accordance with the Clauser formulation. 

Numerical experiments to determine the dependence of < p' > on the level of 
fourth-order dissipation showed that a change in fourth-order smoothing from 0.01 
to 0.05 caused a change of less than 1% in sound pressure level. 

The velocity ratio across the SOFIA cavity shear layer will vary with streamwise 
location. Hence, comparison of the variation of spread rate with velocity ratio is 
shown in Fig. 19b. Three velocity ratios are shown for low Mach numbers with about 
20 points maintained across the layer for all cases. Although the computed spreading 
rates are within the bounds of the experimental data [59], the data point at r = i 
falls below the trend because of the limited entrainment afforded by the inviscid side 
wall treatment. 

5.2.2 Two-Dimensional Rectangular Cavity 

The objective of this two-dimensional cavity case was to demonstrate the prediction 
of self-induced cavity resonance. Validation data are provided by comparison against 
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Rossiter's experiment [15]. Sensitivity of the solution to topology, second-order dis- 
sipation, and turbulence model effects were determined. 

The ca\ity geometry' and grid topology are shown in Fig. 20, where the grid has 
been coarsened for clarity. The test conditions were set as: 

Moo = 0.9, Re L = 1.47 x 10 6 , 1 = 8 in. 

Poo = 0.40 kg/m 3 , = 2.9 x 10 4 N/m 2 

The ratio of cavity length by depth ( L/D ) was 2 for this model. The inflow boundary 
was placed 7.5 L upstream of the cavity leading edge, the outflow boundary 4.5 
L downstream of the cavity trailing edge. The inflow and outflow conditions were 
specified as for the free shear layer cases, and an inviscid wall was placed 5 L above 
the cavity. 



Figure 20: 2-d cavity: topology and grids 


Figure 21 depicts instantaneous Mach number and pressure contours obtained 
during the computation in which the time step size was At = 1.97 ps. Inspection 
of the contours across zonal boundaries indicates that the interpolation process is 
well-behaved for this unsteady flow. The Mach number contours show an instant 
of the time-oscillatory shear layer behavior which is prevalent for these rectangu- 
lar geometries. The pressure contours verify the feedback mechanism postulated by 


CHAPTER 5. RESULTS AND DISCUSSION 



Figure 21: 2-d cavity: instantaneous Mach number and pressure [kPa] contours 
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Rossiter [15], and shown graphically in Fig. 22. Briefly, the cycle begins with the 
propagation of a wave from the aft wall of the cavity to the forward face. Wave 
reflection from the forw r ard wall causes the shear layer to bow' outwards, shedding 
vorticity. The deflected shear layer convccts downstream and induces another cycle. 

The origin of this physical model can be linked to the edge-tone phenomena, where 
a thin planar jet interacts with a wedge. The frequencies at which this feedback is 
reinforced is determined by ambient temperature and Mach number and was first 
quantified by Rossiter. Derivation of this model begins with the assumption that the 
frequency of the vortex shedding is equal to the cavity acoustic field, / = = yA 

The vortical and acoustic field relationships are linked by 

m v X v = L 4- 7„A„ 4- Kuoct' 


L = m a \ a + c T t' 


from which 


__ Uqo ( m - 7) 

= L (i + S) 


( 7 ) 


where the phase lag factor, 7 = 7„ = 0.25, and the normalized convection velocity 
of the perturbations, K = 0.66, are empirically determined constants dependent on 
the geometry and ambient conditions. The integer stage number is given by m = 
m a + m v . Use of the a Mach number scaled by the cavity speed of sound, determined 
by the recovery temperature, offers improved correlation with experiment. The model 
described here is idealized: the shear layer perturbations may be manifested as a 
sinuous motion as opposed to discrete roller vortices depicted in Fig. 22. 

The pressure histories along the cavity walls are depicted in Fig. 23 along with 
the comparison of Rossiter’s data to present results in power spectral density (PSD) 
form. The comparison against experiment shows agreement in frequency at the peak 
magnitudes. Magnitudes are higher for the present case by about 2 dB, however 
for this experiment and other numerical work this has been observed as an effect of 
cavity width. In these studies [15, 36], the sound level was found to be inversely 
related to L/W\ the cavity for Rossiter’s experimental work was of L/\V = 2. Also 
shown in Fig. 23 are the first four stages, m,, predicted by Eq. 7 using both his 
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Figure 22: Rossiter’s feedback model 


vortex convection ratio of K = 0.66 and a K = 0.56 inferred from the treated cavity 
computation discussed below. The PSD for this case was computed using 8192 time 
samples, no zero-padding, and a square window. 

Variation of the second-order dissipation, f 2 , from 0.5 to 0.3 caused no discern- 
able change in the pressure histories. Additionally, in order to test a hypothesis of 
a limited domain of unsteadiness, an isolation zone was implemented as shown in 
Fig. 20. The flow outside the zones of interest was frozen, resulting in a decrease 
of the cavity sound pressure levels by 2%. A final comparison between experimental 
data and numerical results is provided in Fig. 24, where the variation of the mean 
and oscillatory pressures along the cavity walls is shown. The < C' > is computed 
following Rossiter’s reduction of experimental pressure histories 



w'here P and R indicate peak and background pressure levels. The computational 



Power spectral density, PSD, dB Pressure, p, kPa 
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Frequency, f, Hz 


Figure 23: 2-d cavity: pressure history and power spectra comparison 
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Figure 24: 2-d cavity: variation of mean and oscillatory pressures 

results were reduced from assuming p^p = p and p ( \p = p. Use of Reynolds averaging, 
specifically fg = fg, gives a result consistent with the experimental reduction 



Given the difference in spatial dimensions and turbulence modelling uncertainties, 
the trends shown in Fig. 24 appear reasonable for a flow of this complexity. 

Finally, to give a qualitative comparison of numerical and experimental [17] re- 
sults, schlieren images are shown in Fig. 25. Despite Reynolds and Mach number 
mismatches and grid coarseness away from the cavity, the observed and computed 
acoustic radiation patterns are similar. 
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Experiment, Karamcheti (Ref. 17) Present 


Figure 25: Comparison of experimental and numerical schlieren images with knife edge 
vertical 

5.2.3 Two-Dimensional Treated Cavity 

The effect of cavity geometry, particularly modification of the shear layer attachment 
region, is known to possess potential quieting capabilities [42]. The Army Airborne 
Optical Adjunct (AOA), shown in Fig. 26, flight tested several passive and active qui- 
eting methods [94]. The purpose of the present numerical simulations is to determine 
if optical quieting methods, particularly aft ramp treatment and lip-blowing, could 
be accurately simulated. Tangential lip-blowing at the upstream edge of the aperture 
may provide quieting by replenishing the mass entrained by the shear layer from the 
cavity. The quieting provided by aft ramp treatment at the shear layer impingement 
region is discussed below. 

The grid cell size was specified at 0.83" in the streamwise direction, chosen so 
that frequencies up to approximately 400 Hz would be resolved without significant 
numerical dissipation effects. A time step of 44/is was fixed so that CFL ss 1 in 
the streamwise direction within the shear layer. The numerical test conditions were 
matched to flight data: 

Moo = 0.77, Re L = 5.00 x 10 6 , 

poo = 0.262 kg/ Tn 3 , = 1.63 x 10 4 N/m 2 


L = 47 in. 
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Figure 26: U.S. Army Airborne Optical Adjunct [94] 


The initial conditions used the assumption of isentropic recovery to obtain the cav- 
ity temperature while maintaining constant pressure across the aperture. Boundary 
conditions were of the one-dimensional characteristic form: constant mass flow, total 
temperature, and flow angle inflow; constant pressure outflow; and an inviscid wall 
6.4L from the cavity. A characteristic inflow condition was also used for the lip- 
blowing boundary, the flow rate computed using flight data and assuming isentropic 
compression of the ram air utilized in the aircraft. The 100% lip-blowing rate case 
corresponded to a m = 0.42(pu)oo. For the discussion below, computed high and 
low lip-blowing rates refer to 100% and 1% of this mass flow rate. The coarsened 
near-field grids are shown in Fig. 27. 

Comparison of flight data with this planar numerical study is justified by the flight 
test effort toward establishing two-dimensional flow across the apertures of the AOA. 
Use of flow cones and a shear layer rake verified, for the most part, the success of this 
effort. Although the cupola which allows for the cavities is of hemicylindrical form, the 
center of rotation is through the center of the cavity volume. Rather than simulate an 
axisymmetric cavity with an erroneous radius of curvature, the cavity was modelled 
w r ith no spanwise variation in the flow. However, there is a dimensional effect in that 
the mass removed from the cavity by the shear layer entrainment process can only 
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Figure 27: 2-d treated cavity: near field grids 

be replenished at the impingement region. This is in contrast to the mass addition 
mechanism present in three-dimensions, which also includes mass replenishment via 
spanwise structures such as streamwise vortices. 

The mechanism by which an aft ramp reduces the cavity feedback was explained 
by Heller and Bliss [24], Assuming two-dimensional incompressible flow, the region 
immediately surrounding the stagnation point can be treated using a streamfunction 
approach: 

# = axy + ]-by 2 

d* 

u = — = ax + by 

oy 



and the resultant field is shown in Fig. 28. 

Physically, this result can be explained from a force balance normal to a streamline 
approaching the stagnation point. About the stagnation point, the velocity gradient 
across the impinging shear layer creates a pressure gradient. However, there is a 
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T = «xy + -lby* ^ 

Figure 28 : Streamfunction near a stagnation region 


counteracting pressure gradient, pq*/r, due to the differing radii of curvature above 
and below the dividing streamline. 

For a rectangular cavity, the extreme of normal impingement of the shear layer 
onto the aft bulkhead causes further deflection into the cavity. Mass ingestion into 
the cavity causes increased pressure, deflecting the shear layer outwards. With the 
shear region now outwardly deflected, mass expulsion from the cavity reduces the 
cavity pressure, inducing another cycle. Therefore, between the extremes of a normal 
or tangential impingement of the shear layer, a balance of forces may be found. Use 
of a ramp instead of a convex surface at the reattachment region prevents shear layer 
perturbations from inducing instabilities of the type seen in rectangular cutouts. The 
length of the ramp must be large enough to accommodate the magnitude of the 
transverse shear layer excursions expected during operation. 

It is hypothesized that the use of a modestly concave surface at the impingement 
point may provide increased quieting. Improved stability may be provided by the 
following reasoning. As the shear layer is perturbed downwards, the streamlines below 
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the impingement point will have a smaller radius of curvature while the streamlines 
above the shear layer will be less curved. The resultant pressure gradient will drive 
the shear layer upwards. Conversely, as the shear layer deflects upwards the steeper 
tangent angle creates a smaller radius of curvature above, forcing the shear layer 
downwards towards the nominal stagnation point location. 

Computed and flight mean Mach number profiles are compared in Fig. 29 for 
two lip-blowing rates. The quantity <f> indicates the angle from the cupola crest at 
which the data was measured. Figure 29 also shows Mach number contours for the 
two lip-blowing rates above each set of profiles. The Mach number contours are 
instantaneous while the profiles were averaged over 2000 time steps. The difference 
between experiment and computational results on the lower edge of the shear layer 
(r < -2") may be due to blockage in the cavity of the aircraft. This blockage was not 
computationally represented due to the lack of a complete geometry description. The 
difference at the upper aft portion of the shear layer (r > 2") appears to be due to 
blockage in the computational model. Overall, the maximum vorticity as a function 
of x-station is in agreement for both cases. 

Comparisons of power spectra at the aft ramp are shown in Figs. 30 and 31. The 
computed spectra can be seen to be quantitatively and even qualitatively different 
from flight data. The computed result lies more than 15 dB below the data, and 
a peak in the low lip- blowing rate spectra is clearly computed, but is not seen in 
the flight data. The power spectra were computed using 4096 points and a square 
window. Figure 32 shows mean and fluctuating quantities along the cavity w'alls, with 
the available data allowing comparison only along the aft ramp. The computed trend 
in < Cp > is in agreement with measured data, but a large discrepancy in magnitude 
is evident. The large spanwise variation in measured < C' > indicates the existence 
of three-dimensional effects or experimental errors. 

It has been noted from experimental evidence [95] that the frequency of large struc- 
tures in shear layers is independent of axial station and occurs at Strouhal number 
of St = Ui ^ U2 = 0.024 ± 0.003, where 6 is the local shear layer momentum thickness 
and / denotes frequency. This phenomena is corroborated by the reduction of other 
researchers’ data [48, 49, 51, 96, 97], whose results range from about St = 0.02 to 


CHAPTER 5. RESULTS AND DISCUSSION 


Mi .75 



Ms .75 



Figure 29: 2-d treated cavity: instantaneous Mach number contours and mean profiles at 
(a) low and (b) high lip-blowing rate 
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0.03 for incompressible shear layers. 

Gortler’s solution can be used to obtain 6 = 0.036y^x for cr 0 = 11.0, which 
compares favourably to the empirically determined correlation [98] of 9 = 0.034 
Using this relationship along with a compressibility correction [99], the computed 
peak in the AOA solution at 340 Hz corresponds to a Strouhal number of 0.032. For 
comparison, the peak at approximately 1800 Hz in the quieted SOFIA case corre- 
sponds to a St = 0.030, as will be shown in section 5.2.7. It is interesting to note 
that by using Rossiter’s formula, Eq. 7, the frequencies obtained for m = 4 and 5 are 
285 and 360 Hz, respectively. From Fig. 34 it can be observed that m v — 4, implying 
that m a = 1. 


Based upon these observations, it is hypothesized that large scale shear layer struc- 
tures are beginning to be resolved. However, the lack of empirical support from the 
flight data pressure pow r er spectra is at odds with this supposition. The comparison 
is further clouded by the reasonable comparison in < p' > for the low lip-blowing rate 
shown in Fig. 33. The discrepancy may be caused by three-dimensional effects, angle 
of attack sensitivity, or geometry simplifications. 

In this author’s opinion, three-dimensional effects are the most plausible expla- 
nation for the discrepancy. Rockwell [100] noted that for sufficiently large Reynolds 
numbers three-dimensionality reduces coherence in the shear layer. This implies that 
assumption of two-dimensionality for small flow oscillations may be suspect. The evo- 
lution of streamwise-oriented vorticity interacting with the primary vortices would act 
to spread peaks in the reattachment ramp pressure spectra. 


As a final note, Fig. 30 also depicts data, the ordinate scaled by Szshhsbi an( j the 

9oo \tunnel 

abcissa by obtained from an AOA wind tunnel test [101]. The data can 

be seen to agree more closely with the computed results than with flight data, and a 


small peak exists where expected according to the above analysis. 


5.2.4 2-D Treated Cavity: Aero-Optical Effects 

Computation of aero-optical parameters requires the use of the unsteady density field. 
Figure 33 shows the computed and experimental profiles of the root mean square of 
the density fluctuations. Levels of < p' > were computed over a time segment of 
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Figure 32: 2-6 treated cavity: mean and oscillatory coefficients of pressu 
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about 90 ms in increments of 0.44 ms. Using the elapsed time for a particle to 
convect across the aperture at the mean shear layer speed as a characteristic time, 
T c = 77 ^, then the optical computation was taken for about nine T c . In Fig. 33, 
7 is the rake angle from horizontal, with the axis of rotation offset from the cupola 
centerline. Determination of the systematic error band on the experimental result is 
discussed below. The low lip-blowing rate result underpredicts the magnitude of the 
peak in , however the peak location is in fair agreement. The computed results for 
the high lip-blowing rate compare poorly to experiment, possibly due to inadequate 
grid resolution or the increased flow complexity of the merging shear layers. This 
type of active control is presently not a design option for SOFIA, therefore further 
effort toward improvement of the high lip-blowing case was not warranted. 



Figure 33: 2-d treated cavity: profiles with (a) low and (b) high lip-blowing 


Further investigation of the low-blowing rate case revealed the presence of large 
convecting structures associated with the shear layer. Figure 34 shows a contour plot 
°f depicting the growth and propagation of these sinuous motions in the shear 
layer. Also depicted in Fig. 34 is a schematic of the optical model, with the initial 
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and final stations of the optical path integration are given by r 0 and rj. The large 
structures, associated with a 0.03uoo vertical velocity component, are the primary 
contributors of the computed density fluctuations of the shear layer. The speed of 
the waves, as determined from Fig. 35, is 0.56woo, below the value of 0.661(00 inferred 
by Rossiter for rectangular cutouts, yet above the 0.51«oo determined analytically by 
Roscoe and Hankey [102]. 

Chew and Christiansen [50] and Tsai and Christiansen [51], utilizing results from 
computation and experiment, deduced that a free shear layer model of a sinusoidal 
phase delay growing in x w r ould produce results similar to those observed. Figure 35 
displays behavior of a similar nature for the aero-window problem modelled here. 

Comparisons of integrated aero-optical quantities, shown in Fig. 36, reveal slight 
overprediction for the low lip-blowing case and, given the < p' > profiles, expected 
underprediction for the high lip-blowing rate case. Also shown in Fig. 36 are the root 
mean square of the optical path difference fluctuations for two additional integration 
paths. The result for the integration path which extends from ro = —8" to rj = 12" 
displays an increment in < OPD' > of about (7 x 10 -7 )" from the 7" path length 
case which originates at r 0 = —3". The path initialized above the shear layer, from 
ro = 4" to rj = 12", shows a small < OPD' >. Finally, the time averaged optical 
path difference, OPD , can be seen to contribute curvature to the w r avefronts as the 
light propagates through the shear layer. The optical clarity of the shear layer was 
determined using a (3 = 2.584 x 10 -4 , matching the value which was used to reduce 
the experimented data. 

The analytic result for the < OPD' >, which goes like x, is found from [103] 



Derivation of the model, which utilizes time-mean quantities to determine < OPD' >, 
is given by Bogdanoff [103]. This analytic result assumes an index-matched shear layer 
with a sinusoidal n profile, n(y) = sin ( ) . The constants, £ = 0.0091 

and Lson.fn « 1-31 = 1.31(0. 18x[m]), are empirical relations. The virtual origin of 
the shear layer is placed at x = 0" for this analysis. 
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Figure 36: 2-d treated cavity: streamwise variation of optical path difference quantities 
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Reduction of Experimental Data 


The reduction of the data obtained from experiment [104] is noted here to delineate 
the approximations used and estimate systematic error bounds in the optical path 
distortion levels. Values of p' are computed from assumptions of quasi-steady flow: 

h T = c p T (l + ^-A/ 2 ) 
differentiation with respect to t 


RTj 


PP* -Pff , ,7-1 

t 1- uu 

P 2 7 


using (pu)' = pu' + p'u then 


T> 
1 T 

T 


- + (7 - 1)A/ 2 
P 


(puY 

pu 


— [( 7 — 1)A/ 2 + l] — 
L 1 P 


The experimental observations against which the computed results are compared 
assume simultaneously small fluctuations in pressure and total temperature [105, 106], 
resulting in 


i-(— 

P \(7 — 1)A/ ) W 


( 8 ) 


Mean Mach number and density profiles are determined from isentropic relations, 
while teg- is proportional to the voltage fluctuation, f-, obtained from hot film probes. 
The optical path disturbance is then found from [43] 

(OPD'f = 2 (—) P </>’>* l,dr IS) 

V PST P / Jo ' 


where l f is the turbulent eddy size relative to the shear layer width, determined from 
cross correlation data to be typically about 15%. 

The few available independent measurements [106, 107] indicate that pressure 
fluctuations of about 2% of freestream static pressure occur in the shear layer spanning 
the aperture of a quieted cavity geometry. In fact, Hahn [94] reported pressure 
fluctuations of 8% from shear layer rake measurements, however these include the 
dynamic pressure component normal to the orifice as well. Pressure fluctuation levels 
can also be inferred from sound pressure levels in the cavity, observed to be at least 
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130 dB for the AOA case. Shear layer total temperature fluctuations of about 1% 
have also been reported for this Mach regime [107], The present low lip-blowing 
computation found a < p' >« 1% and a < T t >« 0.8% in the shear layer. The 
assumption of < T' t >, < p' >ss 0 in a shear layer is therefore questionable, and is 
used to estimate systematic experimental error bounds. 

The determination of the error in < p' > due to background noise levels begins 
by assuming the passage of a compression wave parallel to the static pressure port in 
the wake rake. Normal reflection of the wave would impart a larger deviation from 
< P' > as computed by Eq. 8. Utilizing Gortler’s free shear layer solution to provide 
u(r), assuming a cavity temperature recovery factor of unity, and holding mean static 
pressure constant through the layer, then p(r) is defined. The sensitivities of ^ to 
f and ^ ^ ± (^ 7 ^ 73+1 and respectively. Using a compression or 

rarefaction wave of strength < p' through the shear layer, then local values of p' 
due to wave passage are defined. This value of p' provides the error bound about the 
value obtained from Eq. 8, which assumes negligible < p' > and < T' r >. Taking 
shear layer pressure fluctuation levels corresponding to 135 dB and a velocity ratio 
T — 0.1, then the systematic error in the density fluctuations is 0.13% at the shear 
layer center. Figure 33 shows the resultant systematic error bars in < p' >. 

From Eq. 9 the value of < p' > is linearly proportional to < OPD' >. The error 
in < OP D' > can be found by using a conservative within-system error of 0.05% in 
gleaned from Fig. 33, plus the systematic error from the above analysis. The 
resultant error bar is plotted in Fig. 36. 

5.2.5 Clean Configuration 

The SOFIA configurations were initialized using the steady solution about a clean, or 
without cavity, geometry. In order to provide a measure of validation, the geometry 
and flow conditions were chosen to replicate the wind tunnel tests: 

Moo = 0.85, Re L = 4.2 x 10 6 , L = 12.6 in. 

Poo = 0.84 kg/m 3 , Poo = 7.7 x 10 4 JV/m 2 , a = 2.5° 
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However, subsequent correction of the wind tunnel data has resulted in the Rei — 
4.0 x 10 6 , about 5% lower than above. The temperature difference from this correction 
causes the sound speed in these computations to be approximately 6% higher than 
experiment. Numerical results obtained for this 7% scale model are discussed below. 

The wind tunnel model without cavity was simulated in order to assess angle 
of attack errors owing to wind tunnel wall effects. Figure 37 compares the present 
pressure coefficient profiles along the crest, side, and bottom of the model with flight 
data [108] and wind tunnel results [16]. Experimental results are shown for both 
untripped and tripped cases; the latter case was used for all subsequent wind tunnel 
testing. The computations specified turbulent walls at all no-slip boundaries. Al- 
though this comparison indicates that the influence of the tunnel wall was small near 
the cavity, pressures along the bottom of the model are shifted, possibly due to the 
effect of the lower wall. A four-order drop in magnitude of 6p\ max was attained for 
this steady case in 2000 steps, using approximately four Cray Y-MP CPU hours. 

5.2.6 Configuration 25 

The geometry shown in Fig. 38 was the initial cavity configuration tested in the wind 
tunnel. This simulation was implemented in order to demonstrate the capture of 
self-excited cavity resonance in three dimensions. The flow conditions were the same 
as used above, the flowfield was initialized from the steady clean case. The stability- 
limited time step used was At = 3.53 ps. This interval size corresponds to a CFL ~ 1 
in the streamwise direction within the shear layer, and a CFL \ max ~ 500. 

Instantaneous Mach number contours in Fig. 39 show the flapping of the shear 
layer and interpolation treatment. Sample pressure histories on the cavity w r alls and 
a comparison of the PSD resulting from the wind tunnel and numerical efforts are 
shown in Fig. 40. The PSD was obtained using 2048 points, a Hanning window, and 
no zero-padding. The predicted frequencies of the dominant tones appear reasonable, 
and the computed dominant tone is within 3 dB of experiment. The magnitudes of 
the computed higher modes are much lower than observed experimentally. 

Estimation of the grid resolution required to maintain a propagating wave of a 
specific magnitude can be deduced from the rectangular two-dimensional cavity and 
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Figure 37: Clean configuration: pressure coefficient comparison 
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Figure 40: Configuration 25: pressure histories and power spectra 
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configuration 25 results. First, wavelength can be estimated by assuming the wave to 
be harmonic at a given frequency and travelling at the local speed of sound. Next, it 
is noted that frequencies around 2 kHz were resolved well in the two-dimensional case, 
in which the grid resolution was such that about 40 points supported the wave. From 
the configuration 25 results, it is seen that only the 700 Hz peak is well resolved, which 
again gives approximately 40 points across the wave for this coarser grid. Although 
numerical damping of the higher frequencies can be expected, most of the energy is 
contained in the lowest frequency mode, as can be seen in the sound pressure level, 




Figure 41: Comparison of sound pressure levels 
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or SPL, comparison of Fig. 41 where 

SPL [dB] = 20 log, 0 

Experimentally observed and computed levels of < p' > in the time domain were used 
in Fig. 41. The SPL for the resonating and quieted geometries obtained numerically 
are in reasonable agreement with experiment. 

5.2.7 Configuration 100 

In 1990, an investigation of SOFIA cavity quieting treatments was performed in the 
NASA Ames 14' x 14' wind tunnel [16]. Of the many geometries tested, configuration 
100 resulted in the lowest sound production levels. This simulation was implemented 
in order to determine if the same level of quieting could be predicted numerically as 
was observed experimentally. As commented on earlier, previous investigations [41] of 
cavity noise suppression have shown aft ramp treatments to be effective by allowing a 
stable shear layer reattachment site. For the SOFIA experiment, this type of geometry 
treatment was found to be quieter than the untreated configuration 25 case by over 30 
dB. Figure 41 summarizes that the proper trends were computed. The flow conditions 
were again initialized from the clean case, and integrated using a stability-limited time 
step size of At = 7.0 6fis. The frequency domain analysis was obtained using 4096 
points, a Hanning window, and no zero-padding. 

For reference purposes, Fig. 42 shows the position and orientation of the telescope 
assembly in the aircraft and the associated coarsened grids. Figure 43 shows the 
topology used in the cavity region, where the grids have again been coarsened for 
clarity. The choice of topology, driven by grid quality and turbulence modelling 
considerations, is similar to those used in the two-dimensional studies. 

Quantitative comparisons were made for this passively quieted cavity geometry in 
terms of shear layer profiles and pressure spectra in the cavity. Since errors in shear 
layer mass entrainment rate would adversely affect the cavity velocity field and hence 
the mean telescope loads, an important validation parameter is the shear layer spread 
rate. Figure 44 depicts mean experimental and computational shear layer Mach 
number profiles. The vertical scale of the profiles is twice that shown in the contour 
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Figure 42: Telescope location and grids 

plot for clarity. Note that the experimental profiles were obtained using a rake, which 
is sensitive only to u, the x-component of velocity. The discrepancy between |u|/c and 
|V|/c was found to be approximately 0.05 at the lower tail of the profile. Figure 44 
indicates reasonable agreement for growth rates, though the profile shapes become 
somewhat different as the shear layer approaches the ramp. This discrepancy may be 
in part due to probe position uncertainty and geometry modifications to allow for the 
probe mechanism. These modifications included removal of the telescope assembly 
and cutting a streamwise slot in the ramp. The difference in spread rates may also 
be due to specification of an overly-large value of Oo- Time averaging of velocities 
was performed over 1000 time steps and the profiles were insensitive to the duration 
of the time-segment used. 

Some measure of qualitative agreement may be gleaned from the instantaneous 
streamlines depicted in Fig. 45, which show a strong cross flow component at the 
aft ramp for this aperture elevation angle. Although oil flow visualization was not 
performed on the configuration 100 experimental runs, a similar aft molding shown 
in Fig. 46 also displays strong cross flow behavior. 

Assessment of the oscillating telescope assembly loads requires the accurate res- 
olution of the unsteady pressure field in the cavity. Comparison of the computed 
and observed spectra at specific locations provides a measure of confidence for the 
computed telescope loads. Toward the estimation of loads, Fig. 47 shows the pres- 
sure history and resultant PSD on the cavity walls. Although the peak levels are in 
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Figure 43: Configuration 100: cavity region topology 
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Figure 44. Configuration 100; instantaneous Mach number contours and mean profiles 



Figure 45; Configuration 100: streamlines along the treated aperture 
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Figure 46: Experimentally [16] observed surface flow pattern 


agreement, the computed spectra can again be seen to drop more rapidly with fre- 
quency than the experimental results. The pressures on the primary mirror, shown 
in Fig. 48, show lower high frequency content with the magnitude of the peak at 1800 
Hz not well resolved. Figure 49 shows a low frequency component at the downstream 
secondary mirror location which was not found experimentally. The discrepancy is 
manifested as the difference between the computed and measured SPL seen in Fig. 41 
at probe 9. 

Scaling to Flight 

The above computations were performed at wind tunnel geometric scale in order to 
allow close comparison with the SOFIA experiments. An early misunderstanding 
resulted in a mismatch of the ambient conditions between the wind tunnel tests and 
the computations. Scaling of optical effects as well as the acoustic frequencies and 
magnitudes from tunnel and computation to flight are obviously of design interest, 
and hence are noted here for completeness and clarity. 
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Figure 47: Configuration 100: pressure history and power spectra on cavity wall 



Power spectral density, PSD, dB Coefficient of Pressure, C 
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Figure 48: Configuration 100: pressure history and power spectra on primary mirror 
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Briefly, the frequency is a function of the cavity length and the recovery temper- 
ature, which affects the acoustic speed in the cavity. The magnitude of the acoustic 
oscillations scale with dynamic pressure, while the optical distortion is proportional to 
density and geometric scale. Model and flight conditions are summarized in Table 2. 
Scaling from wind tunnel to flight results in a decrease in frequency by a multiplica- 



Magnitude 

Quantity 

Units 

Wind tunnel 

Computation 

Flight 

Moc 

- 

0.85 

0.85 

0.85 

Rei 

“ 

4 x 10 6 

4.2 x 10 6 

2.3 x 10 7 

T 

1 oc 

0 K 

286 

322 

217 

Poo 

H7H 

IIktuHI 

0.77 

0.84 

0.289 

L 

m 

0.32 

0.32 

4.6 


Table 2: Wind tunnel, computed, and flight conditions 


tion factor of 16.4, a reduction of the magnitude of pressure oscillations by a factor 
of 11 dB, and an increase in optical distortion by 5.4 times. In contrast, scaling from 
computation to flight results in a frequency reduction of 17.4 times, a decrease in 
fluctuating pressure magnitude of 13 dB, and an increase in optical distortion of 4.9 
times. Generally, the levels of uncertainty in measured quantities is greater than the 
difference between computed and wind tunnel results. 

5.2.8 Configuration 100: Aero-Optical Effects 

The optics code was applied to the computed density field obtained for configuration 
100 from t = 0 to 7.8 ms in the manner depicted in Fig. 50. Ten rays were prop- 
agated through 110 instantaneous density fields in time intervals of A t = 70 .6fis. 
Using the elapsed time for a shear layer structure to convect across the aperture as a 
characteristic time, T c = then the optical measurement was taken for about five 

T c . The results presented here are for a computational plane at approximately the 
cross flow center of the aperture, partly shown in Fig. 50, which will provide only a 
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Figure 50: Configuration 100: instantaneous density field and optical refraction model 

streamwise variation in optical properties. The numerical results are presented com- 
pared to previous analysis [103] and experiment [16] in which shear layer aerodynamic 
measurements were used to infer distortion. 

The levels of fluctuating density were severely underpredicted as compared to 
experiment, as can be seen in Fig. 51. Although peaks in the density fluctuations 
were computed, the highly-ordered shear layer structures similar to those found in 
the AOA study were not observed. Differences may be attributable to grid coarseness 
or within-system errors in measurements, most likely the former. 

The optical wavefront distortion through the configuration 100 aero-window is 
summarized in Fig. 52. The uppermost plot of Fig. 52 shows that the distortion model 
applied through the shear layer alone underpredicts the data determined analytically 
and experimentally. However, the computed trend is generally consistent with the 
data. At the streamwise center of the aperture, the < OPD' > at two additional 
spanwise locations are shown. These points provide an estimate of the crossflow 
variation in distortion levels. 

The center plot of Fig. 52 depicts computed < OPD' > for ray propagation 
originating below the secondary mirror, r 0 = -3.7", and above the shear layer, r 0 = 
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Figure 51: Configuration 100: density fluctuation at crossflow center of aperture 

2.3" . Comparison of the computed results show an increment in < OPD' > below 
the secondary mirror. This distortion increment appears to be caused by a jet of 
re-entrant fluid originating from the shear layer impingement upon the aft ramp. 
Finally, the last plot of Fig. 52 shows that curvature is imparted to the mean optical 
field. The dip in the fluctuating and mean OPD levels at x = 42" is caused by the 
presence of the secondary mirror, in which the index of refraction, n, was fixed at 
unity. 
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Figure 52: Configuration 100: comparison of wavefront distortion 




Chapter 6 
Conclusions 


The objective of this effort was to develop and assess computational methods as ap- 
plied to unsteady perfect gas flows. The developed technology was demonstrated 
through application to two classes of problems in which unsteady effects play a dom- 
inant role. 

6.1.1 Blast- Wave Problem 

The application of two upwind schemes to unsteady, multidimensional problems 
within a structured finite-volume framework has been demonstrated on the viscous 
three-dimensional blast- wave problem. The use of time-conservative differencing and 
an approximate Riemann solver coupled with total variation diminishing methods has 
resulted in time accurate nonoscillatory flowfield resolution. Newton subiterations are 
utilized to reduce the numerical approximations made, such as factorization error and 
the inclusion of only the first-order terms in the formation of the inviscid Jacobian. 
In addition, analysis and application of two flux evaluation methods produced only 
small differences. Finally, for the blast- wave/target interaction problem the effect of 
viscosity was increasingly significant at later times. 

Further efforts to increase the accuracy and efficiency of these methods may be 
directed along the use of nonfactored schemes or implementation on parallel ma- 
chines. Geometries of realistic complexity will require a zonal approach, necessarily 
conservative because of the strongly unsteady compressible flow regimes considered 


105 



CHAPTER 6. CONCLUSIONS 


106 


here. Efficient adaptive grid techniques will reduce the memory and time expense. 
Synthesis of dynamical, structural, and fluid flow effects may provide the capability 
for an interdisciplinary simulation of the physical processes involved in this class of 
problems. 

6.1.2 Cavity Flow Problem 

The work presented here is the initial effort towards development of a cavity flow 
design and analysis tool, specifically tailored for use throughout the SOFIA project 
life. Thus far, this investigation has demonstrated that self-induced cavity resonance 
can be accurately captured for complex geometries modelled using an overset mesh 
topology. Shear layer profiles and resonant behavior are consistent with previous 
analytic and experimental work. Generally, sound pressure levels agree to within 4%. 
Topology treatment has allowed the simple specification of turbulent wall and shear 
layer regions as well as providing a means of isolating the unsteady flow region. 

Improvements in the energy distribution in frequency may be attained by use of 
higher-order spatial approximations or more simply by grid refinement. The use of 
higher order turbulence models should also be investigated. 

6.1.3 Aero-Optical Effort 

Comparison of computed and experimentally observed optical distortion levels showed 
similar trends, albeit with a discrepancy in magnitude. Large structures in the shear 
layer of the two-dimensional quieted cavity resulted in a 25% overprediction of wave- 
front distortion. The three-dimensional results underpredicted phase distortion by 
approximately 50% as compared to experiment. 

Further investigation is required to determine if improved wavefront distortion 
results can be achieved. Quantification of the effects of shear layer flow resolution on 
the optical field is certainly warranted. Improvements in the optical modelling could 
include an empirical model to account for scattering owing to subgrid turbulent scales. 
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6.1.4 Future Directions 

The development of these types of tools partially fulfills the objective of augmentation 
of experimental test programs, possibly eliminating the need to test certain specific 
configurations altogether. The use of the unsteady flowfield information appears fea- 
sible for analysis of the blast-wave/target interaction and cavity flow problem classes. 
However, current limits in computational speed prevent rapid solution throughput. 
Thus, the goal of nesting full unsteady Navier-Stokes methods within the design cycle 
awaits the advent of computers of increased power. 
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Numerical solutions for the tlow Reids about several cav- 
ity conjurations have been computed using the Revnolds- 
averaged Navier-Stokes equations. Comparisons between 
numerical and experimental results arc made* in two- 
dimensions for free shear layers and a rectangular cav- 
ity, and in three-dimensions for the transonic aero-window 
problem of the Stratospheric Observatory For Infrared As* 
'ronomy SOFIA ). Resuits show * hat dominant acoustic 
frequencies and magnitudes ot 'he soil -excited resonant 
cavitv hows compare well with experiment. In addition, 
solution sensitivity to artificial dissipation and grid reso- 
lution levels are determined, furthermore, ootucai path 
distortion due r o the how held :s modelled geometrically 
and is found to match ’lie experimceUai 'rend. 

Introduction 

L ire design of a safe and effective airoorne observatory 
provides motivation tor this study of self-excited open cavi- 
ties. The SOFIA airborne observatory will be a three meter 
class Cassegrain telescope which utilizes a Boeing T o P 
as an observation platform. An artist's concept of the ob- 
servatory, which is currently a cooperative eifort between 
die I’nited Stales’ NASA and Germany's DARA. is shown 
in Fig. 1. Risk of injurv to the crew or damage to the 
platform is present lv oeing reduced dirough experimental 
and computational huid dynamical iCFD) analyses. This 
paper summarizes the analysis methodology developed for 
use in the SOFIA design, of which a detailed outline can 
be found in Atwood and Van Dalsem (1992). 



Fig. 1: Artist s concept of the SOFIA configuration 

Altnougri extensive experimental and numerical studies 
have rwen conducted on the driven cavity problem, ttie ob- 
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jective of this eifort. is to assess muz nemracy ■. -t a g'UaT.d 
and efficient numerical method. I’owarus this goai. the 
following see* ions -'oinpare numerical volutions against an- 
alytical and experimental results for two-dimensionai snear 
! avers and a rcctanguiar cavil y flow Held. Additional vali- 
dation is provided l>v comparison <. f computed SOFIA Row* 
with the experimental data oi Rose and 1 ooley i iTk' • . 

Numerical Method 

General solution of these cavity problems were computed 
using models for the tiuid lieid. die effect ut t erbuience, aim 
die onticai uistortion. the fluid lieid was computed via tiie 
Navier-Stokes equations using die diagonal scheme of Pui- 
liam and Ghaussee >; IDSL, implemented in me overset grid 
framework of Bench. Bulling, and Sieger LTSb u I e<p na- 
tions were integrated through Fuier implicit tune inarm- 
ing and 'econo-order -padai dilferencmg with viscous im- 
permeable waii conditions qjeemeu as no-dip, zero nor- 
mal oressure gradient, and adiaoat ic. Imormatuoii *rans- 
;er across overse* mesn boundaries was implemented us- 
ing triiincar interpolation of die Tuvmd'uit ’.ariabse vector. 
Q — /?, on. pi ' . me. c‘ , 

The eifect oi t urnuictice vns modelled using Me metr.ou 
of Baldwin ami Fornax FG' mouineo 1 id: i variable 
FAar * u t o 1 1 and a .''near !■ ver mocioi. . ue *'d<iy viscos- 
itv in * ne shear .aver was commuted uuu g - it' 1 outer mg ion 
purdori ■ i die Baidwm-iFomax mode], vita ;j ’t “uuai 
to lie) according *o experimental, diougii unite arianm. 
values as compiled by Birch nu Lagers 1 RT JF 

Finaliv, since the speed of ligiit dirougn a gas is a tunc* 
don of the densitv lieid. .mage distortion can * e computed 
using the densitv iiistory ot t. lie snear layer. Ih;s ejtort used 
geometrical ootics m ' Miiunction with tne Loadstone-Dale 
relationship r o ■mmpu r e he mictuating optica; pat!: uit- 
ierence. ■ A i.) T or i . iron which image distortion ’.as 

determined. 


Results 


Free Shear Laver 

Numerical experiments ’-ere performed :>ing a “ wo- 
dimensional mear ovo: *o determine -'‘us:: ivit.-w : mean 
and ' ime- varying e-antith-s 'o nanges m »»•> -w 

fourth-order oissip.t :-.m .cw. -no .::u i nemencui’ 
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Velocity ratio, u.u. 




Fig. 2: (a) Velocity profiles of differing grid resolution com- 
pared to similarity and (b) Variation of spread rate with ve- 
locity ratio 


The results for three grid refinement levels are also shown 
in Fig. 2a along with ' border's similarity solution, where 
the velocity profiles wore taken more than a thousand mo- 
mentum thicknesses downstream of the trailing edge ol the 
plate. L hc solution can he seen to become grid dependent 
when fewer Tuan than 20 points span the layer. Lady vis- 
cosity was observed to grow linearly ;n accordance with 
the Cluuser formulation, while the solution 'vas insensi- 
tive to fourth-order dissipation levels within f c.e range P.01 
to 0.0-5. The Mach ratio for this case was Td/O.fj wir.n 
7 — 20.7, where spreading parameter 7 originates irom tiie 
similarity solution snown in Fig. 2a. I lie spread rate is 
inverseiy related *o c , denoted -Jo w hen tne velocity ot one 
of Tie streams is zero. 

Figure 2h compares ‘he variation >A spread rale with 
velocity ratio, and demonstrates that computed spreading 
rates are 'within Tne hounds of Tie experimental data. 


Two-Dimensional Cavity 

The objective of The two-dimensional cavity »-ompula- 
Ton was to demonstrate the prediction of seif-induced cav- 
ity resonance. Validation data is provided by comparison 
against Rossirer’:? i 1064) experiment. Sensitivity oi ‘he 
solution to tupoiogy, second-order dissipation, and turbu- 
lence model eiFec’s 'were determined. 

Tire rest >'onditions. specified to match experiment, were 
= 6.0. Ret = L .47 x 10°, L = S in., and a cavity length 
by depth. !.■ D, of 2. Characteristic inflow mu outflow 
conditions w>»re 'pecified aiong with a step >ize ut = 

Inspection M Tie computed cavity pressure mstory. 
shown in Fie, 2a. '’oiilirnu the idealized Redback mech- 
anism identified bv Rossi ter. 13 r:eil y. Tie wc le negms with 
Tie propagation oi an acoustic wave trom tne att wail of Tie 
cavity ro Tie forward bulkhead. Wave reflection from the 
forward wall causes the shear layer to bow outwards, shed- 
ding; vorticitv. The deflected shear layer «*onvects down- 
stream and induces another cycle. I he treouencies at which 
Tiis feedback is reinforced is determined by ambient tem- 
perature and Mach number. 

In Tne frequency domain, comparison of Rossi ter s uata 
v> present results indicate agreement in frequency at ’tie 
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Frequency, f, Hz 

Fig. 3: (a) Pressure histories and (b) Power spectra of rect- 
angular cavity 

leak magnitudes, as shown ::i Fig. 2b. Magnitudes arc 1 
higher for f he present case bv about 2 :B. which can ne 
explained from dimensionality argument. 1 lie solution 
was also fowtia to be insensitive to second-order dissipation 
levels within Tie range 0.2 ‘o 0.5. 

.\ umericadv and experimentally determined mean and 
root mean square pressures along Tie '\avity walls are mown 
in Fig. F Fhe ‘rends appear reasonable given Tie flow 
compiexitv. Tie dilference i:i -panai dimensions, and tur- 
bulence modelling assumptions. Finally, the vertical knife 
edge -uTnieren images of Fig. 7 show Tie qualitative agree- 
ment between computed and 'deserved radiation patterns. 
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Fig. -1; Variation of oscillatory and mean pressures 
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b ■ : i n i / ] ? 1 ' i i ■ * : o - : ot t-.u n urn nt ’mu* 

'mqm'ricms. '.'/it h ■ he computed dominant. ‘one wirhin 4 
j 13 ■ '[ experiment. However. -iie masmitudes oi die iiiguer 
modes were much lower diau .ound -experiment ally due ’ o 
grid coarseness. Analysis of the spectra from the resonat- 
ing cavity cases lias shown dial 40 points are required :o 
accurately propagate a travelling wave in the cavity. 
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Fig. 5: Ccnugurat.cn 75: p) pressure history and { b) power 
spectra 

Aldiomm numerical damming ni higiier frequencies can 
,e expected, 'he bulb ■ the energy is contained i:i * he low 
freouent; v modes. os shown in Ftg. 7. [he computed sound 
pressure i-’veis for ‘he resonating and quieted geometries 
ire mas- liable agreement with experiment. 


Configuration 100 

I he ;;ow .d mu' ''nrwiimradon 100 was impiemerred in 
ruer * o uetermme it r lie 'Arne ievei oi 'minting wouid no 
ureuicied numerically as //as found expenmeutaiy/. drevi- 
■ sis investigations *A cavity noise suppression nave mown 
aft ramp treatments to ue eitecdvc ny allowing a table 
-near Fiver reattachinent site. For ’tie oOFIA experiment, 
•ids *voe <c‘ geometrv treatment was louiui ’o i e yiioatr 
dian die um re.atec •' >n;izu radon 11 wise .n«v ver 40 uI3. 
.s mown ;n Fig. A ' sing t.ow c uubjons again initialised 
t an die clean case. die .'taoiiitv- :nite.; ’true - rep one used 
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Fig. 7; Comparison of sound pressure levels 


Figure 7 aiso shows the position ot die bioscope in due 
aircraft and the ait ramp treatment modelled for this case 
aiong with die computed surface how pattern. Figure 
compares the mean experimental and computational snear 
laver Mach profiles, indicating reasonable agreement tor 
growth rates. Die discrepancy in ' tie profile shapes ns * ..e 
shear -aver approaches the ramp may be in parr .;ue to 
an overly high value ^peemed tor prone posh toil en- 
certaintv. or zoo met r’/ moumeations to anuw ;or t no prone 
mechanism. Figure ^ diovs *he general agreement ».i ’tie 
i.ircuu iban< i sued ra. however ‘ ue q>ectra novels wore a siam 
icon to dron more rapiuiy witii frequency tor ’ tie numermai 
as compared to the experimental results. Computed resu.ls 
at die secondary mirror showed a low irequency compormrr 
which was not found experimentally. Fins discrepancy was 
manifested ns a dilicrence in -oimti pressure ’.evms at prove 
location ns shown m Fig. 7. 
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ABSTRACT 

All efficient diagonal scheme implemented in an 
overset mesh framework lias permitted the analysis of 
geometrically complex cavity flows via the Reynolds- 
averaged Navier-Stokes equations. Use of rapid hyper- 
bolic and algebraic grid methods has allowed simple 
specification of critical turbulent regions with an al- 
gebraic turbulence model. Comparisons between nu- 
merical and experimental results are made in two- 
dimensions for the problems of: a backward- facing step, 
a resonating cavity, and two quieted cavity configura- 
tions. In three-dimensions, the flow about three early 
concepts of the Stratospheric Observatory For Infrared 
Astronomy (SOFIA) are compared to wind-tunnel data. 
Shedding frequencies of resolved shear layer structures 
are compared against experiment for the quieted cav- 
ities. The results demonstrate the progress of compu- 
tational assessment of configuration safety and perfor- 
mance. 

NOMENCLATURE 


c 

speed of sound 

Cp 

coefficient of pressure, ^ 

1 

frequency 

K 

ratio of convection by freest ream speed 

( 

mixing length 

L 

characteristic length 

m 

stage number 

in 

mass flow rate 

M 

Mach number 

P 

instantaneous static pressure 

PSD 

power spectral density, dB 


7 

velocity magnitude or dynamic pressure 

Q 

vector of dependent variables 

r 

velocity ratio, liL 

Pc 

St 

Reynolds number 
Strouhal number, — ^ — 
sound power level, dB 

SPL 

t 

time 

u , v , ic 

Cartesian velocity components 


Cartesian physical space coordinates 

O' 

angle of attack 

9 

momentum thickness 

A 

wavelength 

Pt 

eddy viscosity 

P 

density 

(7 

spreading rate parameter 

n 

mean quantity 

< > 

root mean square quantity 

( y 

fluctuating quantity, j — f -F f 

Subscripts 

a 

acoustic 

T 

total quantity 

V 

vortical 

oo 

freest ream quantity 


INTRODUCTION 

The effort reported here describes the progress of 
a computational approach for use in the design of a 
new airborne astronomical observatory. The existing 
Kuiper Airborne Observatory (KAO) has provided two 
decades of unique research capabilities. 1 Figure 1 de- 
picts the proposed successor, SOFIA, which will offer 
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tru times more resolution than the KAO. Assessment 
of the safety ami performance of this large cavity in a 
Hoeing 747SP is the goal of this continuing olfort. To- 
wards this objective of providing design information, 
extensive evaluation of the numerical methods by com- 
parison against experiment is necessary. 

This report describes the status of the validation of 
the computational methods to date. In order to pro- 
vide a somewhat complete overview, both previously 
reported 2 and new results are incorporated herein. The 
two-dimensional cases discussed here include free shear 
layers, 3 a backward-facing step, 4 a resonating cav- 
ity, 0 and two quieted cavities. 0 The computed three- 
dimensional cases are compared to the wind tunnel data 
of Rose and Cooley. 7 The following sections address the 
method used to predict the unsteady (lows, including 
the grid generation and modelling of turbulence. 



Fig. 1: An early conceptualization of SOFIA 


METHOD 

The fluid field was computed via the Navier-Stokes 
equations using the diagonal scheme of Pulliam and 
Chaussee 8 implemented in the overset grid framework 
of Benek, Bulling, and Steger. 9 The equations were 
integrated through Euler implicit time marching and 
second-order spatial differencing with viscous wall con- 
ditions specified as no-slip, zero normal pressure gra- 
dient, and adiabatic. Information transfer across over- 
set mesh boundaries was implemented using trilinear 
interpolation of the dependent variable vector, Q — 
[/?, pu. pi\ pw, e ] 1 . Computations were performed on 
the NAS Cray 2 and Y-MP, the CCF Cray Y-MP, and a 
Silicon Graphics 4D35TG workstation. The flow solver 
cost is 13// s/cell/step on a Cray Y-MP. 


Turbulence Model 

In these computations, the slowly time-varying com- 
ponent of the flow is resolved, while rapid fluctuations 
arc modelled. The algebraic turbulence model of Bald- 


win and Lomax, 10 as implemented by Bulling. 11 is de- 
scribed below using a flow in the (.i\i/) plane. 

The wall-bounded flows use the original Baldwin- 
Loinax model, with the addition of a variable F maz 
cutolL I lie grids are chosen such that a unique wall 
distance is readily available. 

The cavity aperture spanning shear layer uses an 
eddy viscosity developed using F(y) = as sug- 

gested by Baldwin and Lomax 10 for wake regions. This 
results in 


F wake — 


oc 

where specification of C w jt is discussed below and the 
velocity difference is modified to be half the total veloc- 
ity difference between the streams in the specified shear 
layer region 

u, hf = vV 2 + " 2 )„ iai - n /(' (2 + 

The free shear layer model is now given by 

u 2 , f 

IK = P KC' cp C whr f±- ( 1 ) 

[ ^ | rn (i z 

after dropping the KlcbanofF intermit tency function. 
Unmodified model constants K — 0.0 1G8 and C cp = l.G 
are used. 

The magnitude of the eddy viscosity in the free shear 
layer model can be altered with C w k* Estimation of 
the proper value of C u jt begins by using GbrtleEs shear 
layer solution: 
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where u\ and t/ 2 are the velocities of the slow and fast 
streams and ^ is the similarity coordinate, as shown in 
Fig. 2. The spreading parameter cr is inversely related 
to the spreading rate, dbjdr , where b is a measure of the 
shear layer width. The value of the spreading parameter 
when the velocity of one of the streams is zero is <r u . 

GbrtleEs solution can be used to determine the max- 
imum vorticity magnitude as follows: 
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Now, using PrandtTs mixing length assum[)tion and 
scaling laws for jet boundaries, eddy viscosity can also 
be expressed as 
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A' 0 p/;A a 


where AN = tt-t. 

Setting Eqs. 1 and 2 equal results in 
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and only cr 0 remains to be specified. Empirical esti- 
mates of (7 0 range from 9.0 to 13.5/*' 12, 13 For this series 
of cavity flow efforts a q was set to 11.0, resulting in a 
value of C u ,k = 1.91. 


Grid Generation 

Computation of the loads generated by cavity flows 
requires accurate representation of the geometry as well 
as the flowfield. Typically, a significant effort in grid 
generation is required before flow analysis can begin. 
Since matching zone faces arc not required for the over- 
set method used here, recent advances in algebraic 1,) 
and hyperbolic 14 methods can be used. Hyperbolic grid 
generation, which gives good spacing and orthogonality 
control, was used for the wall- bounded regions, while al- 
gebraic grids were used in shear flow regions including 
plumes and wakes. Advantages of this type of grid sys- 
tem include straightforward specification of the turbu- 
lent regions and allowance for independent refinement 
of each zone. This topology also permits the re-use of 
meshes for configuration studies. 


FREE SHEAR LAYER 

Numerical experiments were performed using a two- 
dimensional shear layer to determine sensitivities of 
mean and time-varying quantities to changes in time 
step size, fourth-order dissipation levels, and grid refine- 
ment. In addition, the algebraic turbulent shear layer 
model was partially validated through comparison with 
similarity solutions and experimental data. 

The computational domain for this case is shown to 
scale in Fig. 2a, where inviscid channel side walls and 
characteristic inflow and outflow boundary conditions 
were used. The boundary layers on the splitter plate 



Fig. 2: Shear layer case: (a) velocity profiles of differing 
grid resolution compared to similarity and (b) variation of 
spread rate with velocity parameter 


and the shear la}'er were fully turbulent with a Reynolds 
number based on the mean velocity of the streams and 
the length of the splitter plate given as 6.7 x ID 5 . 

The results for three grid refinement levels are also 
shown in Fig. 2a along with Gortler's similarity solu- 
tion, where the velocity profiles were taken more than 
a thousand momentum thicknesses downstream of the 
trailing edge of the plate. The solution can be seen to 
become grid dependent when fewer than than 20 points 
span the layer. Eddy viscosity was observed to grow lin- 
early in accordance with the Clauser formulation, while 
the solution was insensitive to fourth-order dissipation 
levels within the range 0.01 to 0.05. The Mach ratio for 
this case was 0.2/0. 8 with <j = 20.7. 

Figure 2b compares the variation of spread rate with 
velocity ratio, and demonstrates that computed spread- 
ing rates are within the bounds of the experimental 
data. 

BACKWARD-FACING STEP 

A qualitative estimate of the errors to be expected 
in the three-dimensional flows can be gleaned from the 
backward-facing step problem. The computation was 
compared to the experimental data of Driver and Seeg- 
inillcr, 4 with test conditions matched at 

M vr f = 0.128, Be a = 3.44 x 10* 

p re f = 1.10 k(j/m 3 , Pref =9.11 x 10 4 A"/ III 2 

Characteristic inflow and outflow boundaries were 
used, holding mass flow, total enthalpy, and flow angle 
fixed at the inlet while specifying pressure at the outlet. 
Both til and p ex n were iteratively adjusted to match two 
experimental conditions. First the overall pressure rise 
from .r /// = —4 to .r / II = 30 was matched to experi- 
ment, then m was modified to match the experimental 
momentum thickness at x/H = —4. Finally, viscous 
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adiabatic walls won' specified for doth the top and bot- 
tom of the channel. Tho unusual topology shown in 
Fig. 3 was used to replicate those required by the more 
complex configurations. 



Fig. 3: Backward-facing step: coarsened near field grids 

After dissipation of initial transients, the simulation, 
which was performed using a step size of 15/i.s, was av- 
eraged over 200 steps by increments of 20 steps. Using a 
particle convecting with the shear layer from separation 
to the reattachment point, then 200 stops corresponds 
to three characteristic time increments. The time mean 
field was sampled in /// 10 increments via trilinear in- 
terpolation of a valid donor cell to obtain the profiles 
shown in Fig. 4. 



Mean Mach number contours 



Fig. 4: Backward-facing step: velocity profiles 


tachment occurred at .r/II = G.3 while experimental 
data showed rent t achment at r/H = G . 0 . 

TWO-DIMENSIONAL RESONATING CAVITY 

In order to establish confidence in the numerical 
method, a two-dimensional cavity computation was un- 
dertaken to demonstrate and validate' 1 * self-induced cav- 
ity resonance. 

The test conditions, specified to match experiment, 
were 

M 0 c = 0.9, Rc l = 1.47 x 10°, 1=8 in. 

poo = 0.40 kg/m 3 , p 0 o = 2.9 x 10 1 X/m 2 

and a cavity length by depth, L/D. of 2. Characteristic 
inflow and outflow conditions were specified along with 
a step size of At = 1.07 f is. Power spectra were com- 
puted using 8192 steps following the dissipation of the 
initial transients. 

Inspection of the computed cavity pressure history, 
shown in Fig. 5a, confirms the idealized feedback mccli- 



Fig. 5: 2-d cavity: (a) pressure history and (b) power 
spectra comparison 

anism identified from Rossiter’s experiments. Briefly, 
the cycle begins with the propagation of an acoustic 
wave from the aft wall of the cavity to t lie forward bulk- 
head. Wave reflection from the forward wall causes the 


Agreement in maximum vorticity is seen in Fig. 4, 
however some discrepancies in profile shape exist. Us- 
ing time mean skin-friction levels, the predicted reat- 
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shear layer to how outwards, shedding vortieity. The 
deflected shear layer converts downstream and induces 
another cycle. This coupling of the acoustic and vor- 
tical fields is quantified by Rossiter s empirical model, 
given in Fig. oh, which gives only feedback frequencies. 

In the frequency domain, comparison of Rossiter's 
data to present results indicate agreement in frequency 
at the peak magnitudes, as shown in Fig. 5b. Mag- 
nitudes arc higher for the present case by about 2 clB, 
which can be explained from dimensionality arguments. 
The solution was also found to be insensitive to second- 
order dissipation levels within the range 0.3 to 0.5. 
Figure 5b also shows the resonant modes predicted by 
Rossiter's equation, showing that K = 0.5G gives bet- 
ter prediction of the higher modes. Finally, the vertical 
knife edge schiiercn images of Fig. G show the qualita- 
tive agreement between computed and observed 10 radi- 
ation patterns. 



Experiment, Karamchetl Present 


Fig. 6: Comparison of experimental and numerical 

schlieren images with knife edge vertical 

TWO-DIMENSIONAL TREATED CAVITY 

The effect of cavity geometry, particularly modifica- 
tion of the shear layer attachment region, is known to 
possess potential quieting capabilities. 17 The Army Air- 
borne Optical Adjunct (AOA), shown in Fig. 7, flight 
tested several passive and active quieting methods. 6, 18 
The purpose of the present numerical simulations is to 
determine if optical quieting methods, particularly aft 
ramp treatment and lip-blowing, could be accurately 
simulated. 

The grid cell size was specified as 0.83 in. in the 
streamwise direction, chosen so that frequencies up to 
approximately 400 Hz would be resolved without signifi- 
cant numerical dissipation effects. 19 A time step of 44//s 
was fixed so that CFL % 1 in the streamwise direction 
within the shear layer. The numerical test conditions 
were matched to flight data: 

Moo = 0.77, Rc l = 3.00 x 10 6 

p oo = 0.2G2 kg/m 3 , prx, = 1.63 x 10* N/w 2 

and L = 47 in. The initial conditions used the assump- 
tion of iscntropic recovery to obtain the cavity temper- 
ature while maintaining constant pressure across the 



Fig. 7: U.S. Army Airborne Optical Adjunct 


aperture. A characteristic inflow condition was used for 
the lip-blowing boundary, the flow rate computed using 
flight data and assuming iscntropic compression of the 
ram air utilized in the aircraft. The 100% lip-blowing 
rate case corresponded to a m = 0.42(pu) ::c . For the 
discussion below, computed high and low lip-blowing 
rates refer to 100% and 1% of this mass flow rate. The 
coarsened near-field grids are shown in Fig. 8. 



Fig. 8: 2-d treated cavity: near field grids 

The mechanism by which an aft ramp reduces the 
cavity feedback was explained by Heller and Bliss. 20 
Physically, this result can be explained from a force bal- 
ance normal to a streamline approaching the stagnation 
point. About the stagnation point, the velocity gradi- 
ent across the impinging shear layer creates a pressure 
gradient. However, there is a counteracting pressure 
gradient, pq 2 /r> due to the differing radii of curvature 
above and below the dividing streamline. 

For a rectangular cavity, the extreme of normal im- 
pingement of the shear layer onto the aft bulkhead 
causes further deflection into the cavity. Mass inges- 
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tion into the cavity causes increased pressure, deflect- 
ing th e shear layer outwards. With tin* shear region 
now outwardly deflected, mass expulsion from the cav- 
ity reduces the cavity pressure, inducing another cycle. 
Therefore, between the extremes of a normal or tan- 
gential impingement of the shear layer, a balance of 
forces may be found. Use of a ramp instead of a con- 
vex surface at the reattachment region prevents shear 
layer perturbations from inducing instabilities of the 
type seen in rectangular cutouts. The length of the 
ramp must be large enough to accommodate the mag- 
nitude of the transverse shear layer excursions expected 
during operation. 

Computed and flight mean Mach number profiles are 
compared in Fig. 9 for two lip-blowing rate's. The quan- 


Ms .75 



Mach number 


Fig. 9: 2-d treated cavity: instantaneous Mach num- 
ber contours and mean profiles at (a) low and (b) high 
lip-blowing rate 

tit}' 0 indicates the angle from the cupola crest at which 
the data was measured. Figure 9 also shows Mach num- 
ber contours for the two lip- blowing rates above each 
set of profiles. The Mach number contours arc instan- 
taneous while the profiles were averaged over 2000 time 
steps. The difference between experiment and compu- 
tational results on the lower edge of the shear layer 
(r < — 2") may be due to the computational simplifi- 


cation of an empty cavity. The difference at the upper 
aft portion of the shear layer (r > 2'") appears to be 
due to blockage in the computational model. Overall, 
the maximum vorticitv as a function of x-station is in 
agreement for both cases. 

Comparisons of low lip-blowing rate power spectra 
at the aft ramp are shown in Fig. 10. The computed 



Fig. 10: 2-d treated cavity, low lip-blowing rate: power 
spectra 

result generally falls 20 dB below the flight data, and 
a peak in the low lip-blowing rate spectra at 340 flz is 
clearly computed, but is not seen in the flight data. 
Figure 10 also depicts data, the ordinate scaled by 
and the abcissa by , obtained from 

an AOA wind tunnel test. 21 The data can be seen to 
agree more closely with the computed results than with 
flight data, and a small peak exists at the computed 
frequency peak. The computed high lip-blowing spec- 
tra was broadband, but generally underpredicted flight 
data. 19 Computed power spectra were found using 409G 
points and a square window. 

Experimental data 13 ' 22 ■ 23, 24 indicates that the fre- 
quency of large structures in shear layers occurs at 
Strouhal numbers of St — — ^ — = 0.024 ±0.003, where 
9 is the local shear layer momentum thickness and / de- 


G 





note’s frequency. The computed Strouhal number can 
be estimated using Cottier's solution, given in Fig. 2a, 
to obtain 0 = 0.03Gy-— r for rr n — 11.0, which is com- 
parable to the empirically determined correlation 25 of 
0 = 0.034 jjGr. Specifying r = 0.2 in this relation- 
ship along with a compressibility correction, 2ft the com- 
puted peak in the AOA solution at 340 Hz corresponds 
to a Strouhal number of 0.032. Note that by using 
Rossiter's formula, given in Fig. 5, the frequency ob- 
tained for m = 5 is 3G0 Hz. From Fig. 11 it can be seen 
that four vortical cycles exist (m„ = 4), implying that 
ni a — m — w v = 1. 



Axial station, x, Inches 

Fig. 11: 2-d treated cavity, low lip-blowing rate: instan- 
taneous contours of - £ — 
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The cause of the spike in the power spectrum was 
found from inspection of density fluctuations, which 
revealed large converting shear layer structures. Fig- 
ure 11 shows a contour plot of , depicting the growth 
and propagation of the shear layer structures. The 
large structures, associated with a 0.0 vertical ve- 
locity component, are the primary contributors of the 
computed density fluctuations of the shear layer. The 
speed of t ho structures is 0.5Gu l3O , 19 below the value of 
O.GGmoq inferred by Rossiter for rectangular cutouts, yet 
above the 0.51w.» determined analytically by Roscoe 
and Ilankoy. 28 

Based upon these observations, it appears that large 
scale shear layer structures are being resolved. How- 
ever, the lack of empirical support from the flight data 
pressure power spectra is at odds with this conclusion. 
Three-dimensional effects are a possible explanation for 
t he discrepancy. Rockwell 27 noted that for sufficiently 
large Reynolds numbers three-dimensionality reduces 
coherence in the shear layer. This implies that assump- 
tion of two-dimensionality for small flow oscillations 
may be suspect. In this two-dimensional computation 
the mass removed from the cavity by the shear layer 
entrainment process can only be replenished at the im- 
pingement region. This is in contrast to the mass ad- 
dition mechanism present in three-dimensions, which 
also includes mass replenishment via spanwiso struc- 
tures such as streamwise vortices. The evolution of 


streamwiso-orionted vorticity interacting with the pri- 
mary vortices would act to spread peaks in the reat- 
tachment ramp pressure spectra. 

THREE-DIMENSIONAL RESONATING CAVITY 

An experimental investigation of SOFIA cavity (pil- 
oting treatments was performed in the NASA Ames 
14' x 14' wind tunnel in 1990. 7 The geometry shown 
in Fig. 12 was the initial cavity configuration tested in 
the wind tunnel. This simulation was implemented in 
order to demonstrate the capture of self-excited cav- 
ity resonance in three-dimensions. The simulated flow 
conditions were the same as used in the wind tunnel, 
with the flowfield initialized from a steady clean case. 
A stability-limited time step of At = 3.53 fis was used, 
corresponding to a streamwise CFL « 1 in the shear 
layer, and a CFL\ JUax « 500. 



Fig. 12: Configuration 25: wind tunnel and numerical 
models 

Sample pressure histories on the cavity walls and a 
comparison of the PSD resulting from the wind tunnel 
and numerical efforts are shown in Fig. 13. The PSD 
was obtained using 2048 points and a Hanning win- 
dow to match the treatment of the experimental data. 
The predicted frequencies of the dominant tones appear 
reasonable, and the computed dominant tone is within 
3 dB of experiment. The magnitudes of the computed 
higher modes are much lower than observed experimen- 
tally. 





Fig. 14: Comparison of sound pressure levels 

if the same level of quieting could he predicted numer- 
ically as was observed experimentally. Figure 15 shows 
the topology of the blended aft ramp grids used in this 
simulation of configuration 100. 


Fig. 13: Configuration 25: pressure histories and power 
spectra 

Estimation of the grid resolution required to main- 
tain a propagating wave of a specific magnitude can be 
deduced from the rectangular two-dimensional cavity 
and configuration 25 results. First, wavelength can he 
estimated by assuming the wave to be harmonic at a 
given frequency and travelling at the stagnation speed 
of sound. Next, it is noted that frequencies around 2 
kHz were resolved well in the two-dimensional case, in 
which the grid resolution was such that about 40 points 
supported the wave. From the configuration 25 results, 
it is seen that only the 700 Hz peak is well resolved, 
which again gives approximately 40 points across the 
wave for this coarser grid. Although numerical damp- 
ing of the higher frequencies can be expected, most of 
the energy is contained in t lie lowest frequency mode, 
as can be seen in the sound pressure level, or SPL , 
comparison of Fig. 14. Figure 14 also shows a quieted 
configuration which is discussed below. 


THREE-DIMENSIONAL QUIETED CAVITY 

During the wind tunnel experiments, use of a blended 
aft ramp resulted in the lowest sound production levels. 
This simulation was implemented in order to determine 
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Fig. 15: Configuration 100: cavity region topology 


For the SOFIA experiment, this type of geometry 
treatment was found to be quieter than the untreated 
configuration 25 case by over 30 dB. Figure 14 sum- 
marizes that the proper trends were computed. The 
flow conditions were again initialized from a clean case, 
and integrated using a stability-limited time step size 
of A t = 7.0G//.S. The frequency domain analysis was 
obtained vising 400G points, a Hanning window, and no 
zero-padding. 
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A qualitative comparison of results can be made from 
inspection of the computed instantaneous streamlines 
depicted in Fig. 14 and the oil flow from a similar 
molded geometry, shown in Fig. 1G. A strong cross flow 
component is evident at the aft ramp for t his aperture 
elevation angle, and is seen computationally as well. 



Fig. 16: Experimentally' observed surface flow pattern 

Quantitative comparisons were made for this pas- 
sively quieted cavity geometry in terms of shear layer 
pro files and pressure spectra in the cavity. Since errors 
in shear layer mass entrainment rate would adversely 
affect the cavity velocity field and hence the mean 
telescope loads, an important validation parameter is 
the shear layer spread rate. Figure 17 depicts mean 
experimental and computational shear layer Mach nmn- 
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Mach number 

Fig. 17: Configuration 100: instantaneous Mach number 
contours and mean profiles 

her profiles. The vertical scale of the profiles is twice 
that shown in the contour plot for clarity. Note that the 
experimental profdes were obtained using a rake, which 


is sensitive only to f/. the x-compoiient of velocity. Fig- 
ure 17 indicates reasonable agreement in profile shape, 
however the computed growth rate is lower than seen 
in the data. This discrepancy may be due to geometry 
modifications in the experiment or specification of an 
overly-large value of 

Assessment of the oscillating telescope assembly loads 
requires the accurate resolution of the unsteady pres- 
sure field in the cavity. Comparison of the computed 
and observed spectra at specific locations provides a 
measure of confidence for the computed telescope loads. 
Toward the estimation of loads, Fig. 18 shows the pres- 
sure history and resultant PSD on the cavity walls. Al- 
though the peak levels are in agreement, the computed 
spectra can again be seen to drop more rapidly with fre- 
quency than the experimental results clue to numerical 
dissipation. The spectral peak at 1800 IIz corresponds 
to a Strouhal number of 0.027, using a velocity ratio of 
0 . 1 . 



Fig. 18: Configuration 100: pressure history and power 
spectra on cavity wall 

THREE-DIMENSIONAL QUIETED AFT CAVITY 

Recently a new configuration which places the cavity 
aft of the wing, as shown in Figs. 19 and 20, has been 
proposed as a cost-effective configuration. However the 
cavity response of this aft telescope installation is of 
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Ftg. 19: Aft telescope configuration: modelled geometry 

concern. Limited wind tunnel access lias prevented 
timely resolution of this safety issue through exper- 
imental testing, hence this computation will provide 
an early measure of performance. Flight conditions at. 
41,000 feet and M-^ = 0.85 were used, resulting in a 
Be i = 2.3 x 10'. Powered engines were simulated by 
specification of characteristic conditions for the inlets 
and both fan and core exhaust nozzles. The tail geom- 
etry was not available for this simulation. 
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lty location is shown for reference. Computed pressures 
are generally in agreement with the data, but the com- 
puted gradient aft of the cavity location is more adverse 
than seen in the [light data. This is probably due to dis- 
crepancies in the empennage geometry and the lark of 
tail stabilizers. 
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Fig. 21: Clean configuration: pressure coefficient compar- 


Fig. 20: Aft telescope configuration: cavity grids 

As with the forward cavity computations, a clean air- 
craft flow field was used as an initial condition for the 
cavity computation. Figure 21 compares the computed 
pressure coefficient profiles along the crest, side, and 
bottom of the model with flight data, 29 where the cav- 


Figurc 22 shows the computed PSD on the aft ramp 
as compared to the sealed forward cavity wind tun- 
nel data. The spectra was computed using 1024 points 
stored at 5 At « 0.G8 ms intervals. The computed and 
scaled experimental spectra can be seen to be similar in 
shape and magnitude to about 100 Hz, where numerical 
dissipation becomes significant. The cavity grid cells. 
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Fig. 22: Aft telescope configuration: power spectra 


approximately 3 in. on a side, had been sized so that 
only low frequencies were well resolved. Use of smaller 
cells would capture short wavelength, high frequency 
waves which are unable to excite the massive telescope 
assembly. A Strouhal number of 0.028 is computed for 
the spectral peak seen at 110 Hz. 

Simulation of the distortion of a planar wavetrain 
of light propagating through the cavity aperture 19 re- 
veals the presence of convecting structures, shown in 
Fig. 23. The wavefront distortion, or optical path dif- 
ference (OPD), is proportional to the integrated den- 
sity normal to the aperture. Using f v = j* and Fig. 23, 
it can be verified that the St = 0.028 peak is caused by 
the resolved shear layer structures. The structures con- 
vect with the local downwash velocity field caused by 
the wing. A large spanwise variation of the wavefront 
distortion can also be seen, with maximum < OPD f > 
values at the crossflow edges of the aperture. Since the 
shear layer grid is uniformly distributed in this region, 
this spanwise variation may be indicative of increased 
vortex coherence at the boundaries of the aperture. 

CONCLUSIONS 

The development of an aerodynamic configuration 
analysis method has been summarized. The topology 
used within the overset framework has allowed simple 
and fast mesh generation while maintaining high reso- 
lution in critical flow regions. 

Evaluation of the overset mesh method was per- 
formed through comparison against experiment for sev- 
eral related flows. Sound pressure levels for resonant 
and quieted cavities compared well with experiment. 
However, use of a planar flow assumption for a quieted 
cavity appears to have caused a computed spectral peak 
where flight data indicated broadband behavior. Com- 
puted shear layer structures were found to cause the 



Fig. 23: Aft telescope configuration : wavefront distortion 

spectral peak at a Strouhal number of approximately 
0.03. In three-dimensions, a similar peak was seen in 
both experiment and simulation. Additionally, large 
spanwise density variation of the aperture structures 
was found. 
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Abstract 

The time-varying fluid and optical fields of several 
cavity configurations have been computed on over- 
set mesh systems using the Reynolds-averaged Navier- 
Stokes equations and geometric optics. Comparisons 
between numerical results and Airborne Optical Ad- 
junct (AOA) flight data are made in two-dimensions 
for a quieted cavity geometry with two lip-blowing 
rates. In three-dimensions, two proposed aero-window 
locations for the Stratospheric Observatory For In- 
frared Astronomy (SOFIA) are discussed. The sim- 
ulations indicate that convection of large shear layer 
structures across the aperture cause the blur circle di- 
ameter to be three times the diffraction-limited diam- 
eter in the near-infrared band. 

Nomenclature 

c speed of sound 

dB decibel. 20 log 10 

/ frequency 

h enthalpy 

I intensity 

k wave number, ~L 

A' ratio of convection by freestream speed 

L characteristic length 

m mass flow rate 

M Mach number 

MTF modulation transfer function 

n index of refraction 

OPD optical path difference 

p instantaneous static pressure 

q velocity magnitude or dynamic pressure 

Re Reynolds number 

St Strouhal number, - 1~ 

SR Strehl ratio, -A = exp( — $ 2 ) 

t time 
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T absolute temperature 

T c characteristic time, — 

u,i\w Cartesian velocity components 

x, y<z Cartesian physical space or 

aperture coordinates 
A angstrom, 10” 10 m 

a angle of attack 

3 Gladstone-Dale constant, (n — 1)stp 

0 momentum thickness 

A wavelength 

p density 

a shear layer spreading rate parameter 

$ phase, — ■ ° PD 

( ) mean quantity 

< > root mean square quantity 
( )' fluctuating quantity, / = / -f- /' 
Subscripts 

STP standard temperature and pressure 
T total quantify 

oc freestream quantity 

Introduction 

The study of light propagating through an unsteady 
fluid field has important applications ranging from 
laser weaponry to astronomy platforms. Airborne 
housing of these systems provides mobility, mainte- 
nance, and performance advantages which, in combi- 
nation, can be superior to land or space-based alter- 
natives. However, prediction of the fluid and optical 
behavior of these airborne systems remains a difficult 
problem. 

This report describes the progress of a computa- 
tional approach for use in the design of transonic 
aero-windows. The prediction methodology has been 
driven by the design of the Stratospheric Observatory 
For Infrared Astronomy (SOFIA), the successor to the 
Kuiper Airborne Observatory (KAO), which will offer 
ten times the resolution of the KAO. Figure 1 depicts 
the SOFIA, which will have a 2.5 meter Cassegrain 
telescope mounted in a cavity of a Boeing 747SP. In 
order to numerically assess the safety and performance 
of this platform, extensive evaluation of the computa- 
tional methods by comparison against experiment is 
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Fig. 1: Artist’s concept of the SOFIA configuration 


necessary. 

Many studies of the effect of a fluid field upon 
an optical field have been conducted over the past 
four decades. Many experimental and theoretical ap- 
proaches to the optical distortion problem have been 
investigated; only those of which are pertinent to this 
transonic aero-window problem are summarized here. 
The experimental efforts can be grouped into two cat- 
egories: direct measurement methods and techniques 
based on aerodynamicaJly inferred quantities. Results 
obtained via the latter method are more prevalent be- 
cause of practical difficulties in direct measurement 
techniques. 1 In fact, only aerodynamically inferred dis- 
tortion levels will used for validation of the present 
work. 

Although early experimental and theoretical efforts 
assumed incoherent statistical turbulence, 2 ' 3 4 recent 
studies have begun to examine the effect of shear layer 
structures on electromagnetic field distortion. Using 
a passive scalar field from a direct numerical simu- 
lation, Truman and Lee 0 found an optimum viewing 
angle normal to the hairpin vortices in the homoge- 
neous sheared fluid region. They also found analysis 
via non-refracting geometric optics to be equivalent to 
the parabolized Helmholtz representation of light. Al- 
though this class of studies provides excellent insight 
into the effects of small-scale structure on the electro- 
magnetic field, it is clear that the computational ex- 
pense of such methods precludes their near-term use 
for the problems under consideration here. 

The study of large scale structures in shear layers 
has been an active topic of research since they were 
observed by Brown and Roshko in 1974. 6 Only re- 
cently has the effect of these structures on the optical 
field been studied. In 1990, Chew and Christiansen 7, 8 
experimentally observed the effect of shear layer struc- 
tures on beam propagation. Tsai and Christiansen 9 
used an Euler simulation to determine the optical char- 
acteristics of a perturbed free shear layer. It was hy- 
pothesized that the effect of the vortical structures on 
the optical field could be modelled by a growing sinu- 
soidal phase plate. 

The numerical modelling of the optical effect of a 
cavity-spanning shear layer was presented by Cassady, 
Birch, and Terry 10 in 1987. They found their two- 


dimensional solution to result in poor prediction of 
optical distortion. Farris and Clark 11 12 used time- 
mean quantities and empirical evidence to ascertain 
the fluctuating density levels required for optical anal- 
ysis. 

The present effort attempts to determine what por- 
tion of the optical path distortion can be resolved using 
cell sizes required to obtain an accurate flowfield so- 
lution. Towards this end. computed optical distortion 
levels are compared to flight or wind tunnel measure- 
ments for two- and three-dimensional quieted cavities. 

Previous reports 1314 have described the method 
development for two-dimensional free shear layers, a 
backward-facing step, and a rectangular cavity. 15 Com- 
parison of the computed cavity case with Rossiter's 
data showed agreement in the dominant resonant 
peaks to within 5 dB. The computed and experimen- 
tally observed pressure loading trends were similar 
along the cavity walls. In three-dimensions, rectangu- 
lar and treated quiet cavity solutions were computed 
and compared to experiment. 16 Sound pressure lev- 
els along the cavity bulkheads for both the resonating 
and quieted geometries were found to be in agreement. 
However, although the power spectra of the experi- 
ment and computation were similar at low frequen- 
cies, numerical dissipation caused a rapid decrease in 
energy content at high frequencies. 

Although the optical model has been described pre- 
viously, 13 this paper documents the extension of the 
model and provides new validation information. The 
following sections address the method used to predict 
the unsteady flows and the resultant optical distortion. 
Analysis of the aperture fluid and optical fields for 
AOA and SOFIA configurations are presented. For the 
two-dimensional AOA geometry, time- varying density 
fields and optical path lengths are shown. Short and 
long exposure far-field diffraction patterns are com- 
puted for a three-dimensional aft cavity SOFIA con- 
cept. 


Approach 

Solution of the aircraft and cavity flowfields were 
computed using models for the fluid field, the effect 
of turbulence, and the optical distortion. A diagonal 
scheme 1 ' was used for the solution of the Reynolds- 
Averaged Navier-Stokes equations, implemented in an 
overset grid framework. 18 Euler implicit time integra- 
tion and second-order spatial differencing was used, 
with viscous impermeable wall conditions specified as 
no-slip, zero normal pressure gradient, and adiabatic. 
Information transfer across overset mesh boundaries 
was implemented using trilinear interpolation of the 
dependent variable vector, Q = [p, pu, pi\ pu\e] r . Al- 
gebraic turbulence models were used, implemented 
with a variable F ma x cutoff for wall-bounded flows and 
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a shear layer model for the cavity aperture region. 19 20 
The flow solver cost was 13 // a*/ rell/iteration on a single 
head of the Numerical Aerodynamic Simulator (NAS) 
Cray Y/MP-832. 

Generally, a significant effort in grid generation 
is required before flow and optical analysis can be- 
gin. However, recent advances in algebraic 21 and 
hyperbolic 22 methods have enabled rapid discretiza- 
tion of complex geometries. Hyperbolic grid genera- 
tion, which provides spacing and orthogonality control, 
was used for the wall-bounded regions, while algebraic 
grids were used in shear flow regions including plumes 
and wakes. This choice of topology allows simple spec- 
ification of turbulent regions and also permits the recy- 
cling of meshes, useful for configuration changes such 
as cavity positioning. 

The optical computations documented here use a 
refracting-ray method, reported on earlier, which is 
limited to studying the effects of the resolved large- 
scale structures. 23 The method tesselates a structured 
grid into tetrahedra and uses piecewise mean indices 
of refraction for each of the tetrahedra. Indices of re- 
fraction were computed using n = 1 + {3 — , where 
the Gladstone-Dale constant, /?, can be found using 
the Cauchy formula. 

Assessment of the optical performance of an aero- 
window begins by specification of the ray initialization 
plane. Integration of the optical path length through 
the aero-window is then performed along the ravs at 
specified time increments in both the streamwise, x, 
and crossflow, y, aperture directions. The resultant 
OPD(t,x>y) can be used in a complex aperture func- 
tion of the form P = A Assuming no loss in 

transmittance, then the wave amplitude .4 is unity in 
the aperture and zero elsewhere, while the phase of the 
wave is computed from OPD(t,x,y) = . The au- 

tocorrelation of the complex aperture function, P * P, 
gives the far field diffraction pattern, computed using 
a two-dimensional Fourier transform. Time averag- 
ing successive short-period diffraction patterns gives a 
long exposure result. Integration of the intensity of 
this resultant long or short exposure diffraction pat- 
tern gives the area for a specified encircled energy level. 
From this area the equivalent blur-circle diameter due 
to the resolved fluid scales can be found. Inclusion of 
the effects of small-scale turbulence could be incorpo- 
rated into the computation of long exposure blur cir- 
cles by multiplication of the above modulation transfer 
function (MTF) with a turbulence MTF. 3, 24 

Results and Discussion 

Aero-optical simulations of the U.S. Army Airborne 
Optical Adjunct (AOA) and the SOFIA configurations 
are discussed below. Information pertaining to the 
computation and analysis of the unsteady flowfields. 


including grid resolution and turbulence modelling, is 
given elsewhere. 14 23 


2-D AOA Cavities 

Data available from flight tests of the AOA 25 shown 
in Fig. 2, provides valuable validation information for 



Fig. 2: U.S. Army Airborne Optical Adjunct 


the present simulations. These two-dimensional nu- 
merical simulations were used to determine if optical 
quieting methods, particularly aft ramp treatment and 
lip-blowing, could be accurately simulated. The flow 
about the geometry, depicted in Fig. 3, was computed 
in conjunction with two lip-blowing rates for the for- 
ward aperture only. The 100% lip- blowing rate case 
corresponded to a m = 0.42 (pu)^. For the discussion 
below, computed high and low lip-blowing rates refer 
to 100% and 1% of this mass flow rate. 


M = .75 



Fig. 3: AOA case: instantaneous Mach contours 


Computed and flight mean Mach number profiles 
are compared in Fig. 4 for two lip-blowing rates, with 
instantaneous Mach number contours shown above 
each set of profiles. The quantity <p indicates the angle 
from the cupola crest at w'hich the data was measured. 
Overall, the Mach number contours, which were aver- 
aged over 2000 time steps, show agreement in the max- 
imum vorticity as a function of x-station. However, 
differences exist between experimental and computa- 
tional results at the low-speed edge of the shear layer. 
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Fig. 4: AOA case: instantaneous Mach number contours 

and mean profiles at (a) low and (b) high lip-blowing rate F '*■ 5: A0A case: P ower spectra, low lip-blowing rate 


which may be due to blockage in the cavity of the air- 
craft that was not computationally represented. The 
discrepancy at the upper aft portion of the shear layer 
appears to be due to blockage in the computational 
model. 

Comparisons of pressure spectra at the aft ramp for 
the low lip-blowing rate are shown in Fig. 5. The com- 
puted spectra can be seen to be quantitatively and 
even qualitatively different from flight da* a. The com- 
puted result lies more than 15 dB belo\s ue data, and 
a peak in the low lip-blowing rate spectra is clearly 
computed, but is not seen in the flight data. The high 
lip-blowing rate spectra were similar to that shown in 
Fig. 5, albeit without the spectral peak at 340 Hz. 23 

It has been noted from experimental evidence 26 that 
the frequency of large structures in shear layers is inde- 
pendent of axial station and occurs at Strouhal num- 
ber of St = = 0.024 ± 0.003, where 9 is the 

Uf “Ttij 

local shear layer momentum thickness and / denotes 
frequency. This phenomena is corroborated by the re- 
duction of other researchers' data, 6 7 9 27 • 28 whose 
results range from about St = 0.02 to 0.03 for incom- 
pressible shear layers. 

Momentum thickness can be estimated using 
Gortlers solution, giving 0 = 0.036y^a* for <r 0 = 11.0, 
which compares favorably to the empirically deter- 
mined correlation 29 of 0 = 0.034|^x. Using this re- 
lationship along with a compressibility correction, 30 
the computed peak in the AOA solution at 340 Hz 


corresponds to a Strouhal number of 0.032. For com- 
parison, peaks were found in SOFIA experiments and 
computations at approximately St = 0.028. 

Based upon these experimental observations, it is 
hypothesized that large scale shear layer structures are 
being resolved. However, the lack of empirical sup- 
port from the flight data pressure power spectra is at 
odds with this conclusion. The comparison is further 
clouded by the reasonable comparison in < p r > for the 
low lip-blowing rate shown in Fig. 6, and the presence 
of the organized structures shown in Fig. 7. 
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Fig. 6: AOA case: -f- ^ profiles with (a) low and (b) 
high lip-blowing 
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Three-dimensional efforts are a possible explana- 
tion for the discrepancy between computation and 
flight. Rockwell 31 noted that for sufficiently large 
Reynolds numbers three-dimensionality reduces coher- 
ence in the shear layer. This implies that an error in 
the assumption of planar How for the small flow os- 
cillations considered in these quieted cavity configu- 
rations. The evolution of streamwise-oriented vortic- 
ity interacting with the primary vortices would act to 
spread peaks in the reattachment ramp pressure spec- 
tra. As a final note. Fig. 5 also depicts data, the ordi- 
nate scaled b\ {Qoc ) flight / {Qoc)tunnei and the abcissa 
by (Coc/ L) flight/ (Coc-/ L)tunneh obtained from an AOA 
wind tunnel test. 32 The wind tunnel data shows that 
a small peak exists where expected according to the 
above analysis. 



Fig. 7: AOA case, low lip-blowing rate: instantaneous 
contours of (a) and (b) with schematic of optical 
model 

Using the computed unsteady density field, aero- 
optical effects can be determined. Figure 6 compares 
the computed and experimental profiles of the root 
mean square of the density fluctuations. Levels of 
< p’ > were computed over a time segment of about 
90 ms in increments of 0.44 ms. Using the elapsed 
time for a particle to convect across the aperture at the 
mean shear layer speed as a characteristic time, then 
the optical computation was taken for about nine T c . 
In Fig. 6, 7 is the rake angle from horizontal, with the 
axis of rotation offset from the cupola centerline. De- 
termination of the systematic error band on the exper- 
imental result is discussed below. The low lip-blowing 
rate result underpredicts the magnitude of the peak 
in < ? > , however the peak location is in fair agree- 
ment. The computed results for the high lip-blowing 
rate compare poorly to experiment, possibly due to 
inadequate grid resolution and/or the increased flow 
complexity of the merging shear layers. This type of 
active control is presently not a design option for the 
SOFIA, therefore further effort toward improvement 
of the high lip-blowing case was not warranted. 

Further investigation of the low-blowing rate case 
revealed the presence of large convecting structures as- 


sociated witli the shear layer. Figure 7 shows a contour 
plot of . depicting the growth and propagation of 
these shear layer structures. Also depicted in Fig. 7 is 
a schematic of the optical model, with the initial and 
final stations of the optical path integration given by r 0 
and r f. The large structures, associated with a 0.03^^ 
vertical velocity component, are the primary contribu- 
tors of the computed density fluctuations of the shear 
layer. The speed of the waves, as determined from 
Fig. 8. is O.oCuoc, below the value of 0.66*/^ inferred 
by Rossi ter la for rectangular cutouts, yet above the 
0.5 determined analytically by Roscoe and Han- 
key. 33 
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Fig. 8: AOA case: contours of OPD'(xJ)[ in.] along 
aperture, low lip-blowing rate 

In 1990. Chew and Christiansen 8 and Tsai and 
Christiansen 9 deduced that a free shear layer model 
of a sinusoidal phase delay growing in x would pro- 
duce results similar to those observed in both compu- 
tation and experiment. Figure 8 displays behavior of a 
similar nature for the aero- window problem of concern 
here. 

Comparisons of integrated aero-optical quantities, 
shown in Fig. 9, reveal slight overprediction for the 
low lip-blowing case and, given the < p* > profiles, 
expected underprediction for the high lip-blowing rate 
case. Also shown in Fig. 9 are the < OPD f > for two 
additional integration paths, shown to demonstrate 
that most of the optical distortion is caused by the 
shear layer. 

The result for the integration path which extends 
from r 0 = —8" to r/ = 12" displays an increment in 
< OPD' > of about (7 x 10 -7 )" from the 7 " path 
length case which passes through the shear-layer alone. 
The path initialized above the shear layer, from r 0 = 
4" to rf = 12", shows a small < OPD f > indicating 


o 


that < p > above the shear layer is small. Finally, 
the time mean optical path difference, TTPD , can be 
seen to contribute curvature to the wavefronts as the 
light propagates through the shear layer. The optical 
clarity of the shear layer was determined using a 3 = 
2.584 x 10~ 4 . matching the value which was used to 
reduce the experimental data. 



Fig. 9: AOA case: streamwise variation of optical path 
difference quantities 


The analytic result for the < OPD ( >, which goes 
like j\ is found from 34 



Derivation of the model, which utilizes time-mean 
quantities to determine < OPD' >, is given by 
Bogdanoff. 34 This analytic result assumes an index- 
matched shear layer with a sinusoidal n profile. n(y) = 

sin ( o r 2Ty — V The constants, — = 0.0091 
and L$Q n m a 1.315*. = 1.31(0.18x[m]), are empirical 


relations. The virtual origin of the shear layer is placed 
at jt = 0 to obtain this analytic result. As shown in 
Fig. 9, the analytic gradient in < OPD ! > is in good 
agreement with flight data. 


Reduction of Experimental Data 

The reduction of the data obtained from experi- 
ment 30 is noted here to delineate the approximations 
used and estimate systematic error bounds in the op- 
tical path distortion levels. Values of p are computed 
from assumptions of quasi-steady flow: 


h T = c p T^ l + ^-A/ 2 ) 
Differentiation with respect to t gives 


= 



+ uu' 



Using (pu) f = pu f 4- p'u then 


V 

1 T 

T 


^- p - + - [ ( , - dm' + 1] e - 

5 PU 1 J p 


The experimental observations against which the com- 
puted results are compared assume simultaneously 
small fluctuations in pressure and total tempera- 
ture, 36, 3 ‘ resulting in 


P 


(— 4 =, 

Vh-i W 



(pu)' 

pu 


( 1 ) 


Mean Mach number and density profiles are deter- 
mined from isentropic relations, while is propor- 
tional to the voltage fluctuation, obtained from 
hot film probes. The optical path disturbance is then 
found from 1 


(OPD 1 ) 2 = 2 ( — — — T f L < p' > 2 l r dr (2) 
\PSTPj Jo 

where ^ is the turbulent eddy size relative to the shear 
layer width, determined from cross correlation data to 
be typically about 15%. 

The few available independent measurements 37 38 
indicate that pressure fluctuations of about 2%. of 
freestream static pressure occur in the shear layer 
spanning the aperture of a quieted cavity geometry. 
In fact, Hahn 25 reported pressure fluctuations of 8% 
from shear layer rake measurements, however these in- 
clude the dynamic pressure component normal to the 
orifice as well. Pressure fluctuation levels can also be 
inferred from sound pressure levels in the cavity, ob- 
served to be at least 130 dB for the AOA case. Shear 
layer total temperature fluctuations of about 1% have 
also been reported for this Mach regime 38 The present 
low lip-blowing computation found a < p* >» 1% and 


6 



a < Tj 0.8V: in t he shear layer. The assumption 
of < T'j >. < // 0 in a shear layer is therefore 

questionable, ami is used to estimate systematic ex- 
perimental error bounds. 

The determination of the error in < p l > due to 
background noise levels begins by assuming the pas- 
sage of a compression wave parallel to t he static pres- 
sure port in the wake rake. Normal reflection of the 
wave would impart a larger deviation fr a < p f > 
as computed by Eq. 1. Utilizing Gortler's free shear 
layer solution to provide u(r), assuming a cavity tem- 
perature recovery factor of unity, and holding mean 
static pressure constant through the layer, then p(r) 
is defined. The sensitivities of 4 to i and are 

P p Tt 

± ( 7 -da/ j + i and VTV’ + V- respectively. Using 

a compression or rarefaction wave of strength < p' > 
through the shear layer, then local values of p f due to 
wave passage are defined. This value of p f provides 
the error bound about the value obtained from Eq. 1, 
which assumes negligible < p f > and < V T >. Taking 
shear layer pressure fluctuation levels corresponding 
to 135 dB and a velocity ratio r = 0.1, then the sys- 
tematic error in the density fluctuations is 0.13% at 
the shear layer center. Figure 6 shows the resultant 
systematic error bars in < p' >. 

From Eq. 2 the value of < p’ > is linearly propor- 
tional to < OPD' >. The error in < OPD' > can be 
found by using a conservative within-system error of 
0.05%) in — gleaned from Fig. 6. plus the systematic 
error from tfie above analysis. The resultant error bar 
: plotted in Fig. 9. 

Forward Cavity SOFIA 

Wind tunnel tests, completed in 1990. allowed vali- 
dation of cavity acoustic response and optical charac- 
teristics. The cavity environment results are summa- 
rized in Fig. 10 for both resonant and quieted config- 
urations. 13 Aerodynamic measurements in the shear 
layer were used to infer optical quantities for the qui- 
eted cavity, configuration 100. 

Following computation of the unsteady flow, the op- 
tics code was applied to the computed density field 
obtained for configuration 100 from t = 0 to 7.8 ms. 
Ten rays were propagated through 110 instantaneous 
density fields in time intervals of At = 70.6//S. The op- 
tical measurement was taken for about five T c , again 
using the elapsed time for a shear layer structure to 
convect across the aperture as the characteristic time. 
The forward cavity SOFIA results presented here are 
for a computational plane at approximately the cross 
flow center of the aperture, which will provide only a 
streamwise variation in optical properties. The numer- 
ical results are presented compared to previous anal- 
ysis 34 and experiment 16 in which shear layer aerody- 



Fig. 10: Forward cavity SOFIA cases: Comparison of 
sound pressure levels 

namic measurements were used to infer distortion. 

Comparison of experimental and numerical density 
fluctuations are shown in Fig. 11, where < p f > can be 
seen to be severely underpredicted. Although peaks 
in the density fluctuations were computed, the highly- 
ordered shear layer structures similar to those found m 
the AOA study were not observed. Differences may be 
attributable to grid coarseness or within-system mea- 
surements errors. 
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Fig. 11: Configuration 100: density fluctuation at cross- 
flow center of aperture 

The optical wavefront distortion through the con- 
figuration 100 aero-window is summarized in Fig. 12. 
Figure 12a shows that the distortion model applied 
through the shear layer alone underpredicts the data 
determined analytically and experimentally. However, 
the computed trend is generally consistent with the 
data. At the streamwise center of the aperture, the 
aerodynamically inferred < OPD ' > at two additional 



spanwiso locations arc shown. These points provide an 
estimate of the crossflow variation in experimental dis- 
tortion levels. 
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Fig. 12: Configuration 100: comparison of wavefront 
distortion 

Figure 12b depicts computed < OPD' > for ray 
propagation originating below the secondary mirror, 
ro = —3. 7". and above the shear layer, r 0 = 2.3". 
Comparison of the computed results show an incre- 
ment in < OPD' > below the secondary mirror. This 
distortion increment appears to be caused by a jet of 
re-entrant fluid originating from the shear layer im- 
pinging on the aft ramp. Finally, Fig. 12c shows that 
curvature is imparted to the mean optical field. The 
dip in the fluctuating and mean OPD levels at x = 42" 
is caused by the presence of the secondary mirror, in 
which the index of refraction, n, was fixed at unity. 

Aft Cavity SOFIA 

Forward placement of the telescope in a favorable 
pressure gradient region has an advantage in terms of 
an optically thin boundary layer. However, the fuse- 
lage moldline and structural complexities forward of 


the wing present considerable manufacturing difficul- 
ties. An alternative site for the telescope aft of the 
wing reduces the modification costs and permits the 
use of a larger usable cavity volume. However, an aft 
cavity site has potential problems of scattered light 
emitted from engines and plumes, an optically thick 
boundary-layer, unknown cavity response, and possi- 
bly poor empennage flow behavior at off-design condi- 
tions. 

Figure 13 depicts the simplified geometry used to 
address some of these concerns; horizontal and vertical 
stabilizer geometry was unavailable for this simulation. 
Details of this flight condition simulation are available 
elsewhere. 14 



Fig. 13: Aft SOFIA case: Surface C p and plume tem- 
perature contours 

The acoustic response of the aft cavity is compared 
to scaled data from the forward cavity experiment in 
Fig. 14. The computed result is taken from a location 


130 



Fig. 14: Aft SOFIA case: Power spectra comparison 

on the aft ramp, while the experiment is from a loca- 
tion within the cavity. The agreement of spectra is rea- 
sonable to about 100 Hz, above which grid coarseness 
dissipates energy rapidly. 14 Figure 14 shows peaks at 
60 and 110 Hz, the latter corresponding to a Strouhal 
number of 0.028. 

During the aft cavity SOFIA computation the en- 
tire aperture density field was saved in increments of 
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0.68 m.s\ ever y five flow solution steps. Using this den- 
sity field, propagation of a plane wave through the 
aperture revealed variations in the wavefront distor- 
tion, as shown in Fig. 15. These ordered variations 
in OPD . indicative of shear layer structures in the 
aperture, impart the aft ramp at a frequency of 110 
Hz. giving a St — 0.028. Figure 15 shows maximum 
distortions of about one wavelength. ±Ap, with a re- 
sultant maximum < OPD' > of approximately 0.7A^. 
Computation of the OPD{t<x s y) was performed for a 
64 x 64 array of rays normal to the aperture and ini- 
tialized just above the secondary mirror. The optical 
integration was performed over 8 T c . 



Fig. 15: Aft SOFIA case: Sample wavefront distortion 
history 


Using these phase distortion levels, far field diffrac- 
tion patterns were then computed. Figure 16 depicts 
the diffraction-limited Airy pattern for reference, and 
both instantaneous and time-averaged exposures. The 
instantaneous exposure pattern shows some evidence 
of speckle, with a large reduction in central peak in- 
tensity. This spreading of energy is manifested in the 
computed Strehl ratio of 0.34. 

Finally, Fig. 17 shows that the large scale structures 
in the shear layer cause the equivalent 80% blur cir- 
cle to be three times the diameter of the diffraction- 
limited case. However, as can be seen in Fig. 16. the 
blurring in the streamwise direction is worse than in 
the crossflow direction. Note that the because the 
small scale fluid motion is modelled when using the 
Reynolds-averaged Navier-Stokes equations, the com- 
puted blur circle is much smaller than actually ob- 
served. 39 
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Fig. 16: Aft SOFIA case: Far field diffraction patterns 
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Fig. 17: Aft SOFIA case: Normalized integrated intensity 
distributions 

Conclusions 

Computations of quieted cavity configurations have 
shown convection of large scale flow structures across 
the aperture. The shedding frequency of these struc- 
tures compare reasonably well with experimentally de- 
termined shear layer Strouhal numbers. The com- 
puted results indicate that three-dimensional effects on 
the shear layer spanning a quieted cavity can be signifi- 
cant. The differences in two- and three-dimensional re- 
sults are manifested in the spectra of the pressure loads 
and in the magnitude of the optical wavefront distor- 
tion. Since the primary contributors to the computed 
OPD' were the large scale structures, computations of 
the Strehl ratio were found to be reasonable. However, 
because the small scale fluid motion is modelled, the 
blur circle diameter is significantly underpredicted. 

Further improvements to the prediction of optical 
performance may be found from investigation of shear 
regions with direct Navier-Stokes methods 5 or use of a 
turbulence MTF. 24 Finally, although low Reynolds and 
Mach number number experiments show large struc- 
ture formation, direct OPD measurements at realistic 
conditions would be useful for validation of numerical 
aero-optical studies of this type. 
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Abstract 

Coupling of the Reynolds-averaged Navier-Stokes equations, rigid-body dynamics, and a pitch at- 
titude control law is demonstrated in two- and three-dimensions. The application problem was the 
separation of a canard-controlled store from an open-flow rectangular cavity bay at a freestream Mach 
number of 1.2. The transient flowfield was computed using a diagonal scheme in an overset mesh frame- 
work, with the resultant aerodynamic loads used as the forcing functions in the nonlinear dynamics 
equations. The proportional and rate gyro sensitivities were computed apriori using pole placement 
techniques for the linearized dynamical equations. These fixed gain values were used in the controller 
for the nonlinear simulation. Reasonable comparison between the full and linearized equations for a 
perturbed two-dimensional missile was found. Also in two-dimensions, a controlled store was found 
to possess improved separation characteristics over a canard-fixed store. In three-dimensions, trajec- 
tory comparisons with wind-tunnel data for the canard-fixed case will be made. In addition, it will 
be determined if a canard-controlled store is an effective means of improving cavity store separation 
characteristics. 


Nomenclature 


C m 

coefficient of pitching moment 

DOF 

degree of freedom 

I 

body inertia 

K a 

proportional gain 

K r 

rate gyro sensitivity 

L 

characteristic length 

m 

body mass or stage number 

M 

Mach number 

q 

dynamic pressure 

Q 

vector of dependent variables 

Re 

Reynolds number 

s 

Laplace operator 
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/ time 

T absolute temperature 

u.v.w Cartesian velocity components 

x, //. 2 Cartesian physical space 

ot angle of attack 

£> effector deflection angle 

p density 

( 7 shear layer spreading rate parameter 

r canard servo characteristic time 

6 body attitude 

u) vorticity 

( ) time derivative 

Superscripts 

n time step 

Subscripts 

c commanded, canard 

t tail 

oo freest re am quantity 

Introduction 

The design of high-performance aircraft has in the past typically been compartmentalized by dis- 
cipline, most commonly structures, fluids, and controls. However, as performance demands escalate, 
aircraft systems have become increasingly interrelated. 1, 2 Therefore, there is a need to investigate the 
optimization of the aircraft in its entirety, not simply by evaluation of its sub-systems. 

This type of multidisciplinary analysis is currently accomplished with simplified physical models, 
such as panel flow and modal structural codes. However, in critical regions of the flight envelope these 
linear methods can fail, leading to the necessity for higher-order models. 3, 4 Obviously, this increased 
physical fidelity comes with a high computational price, and hence fewer design cycles are permitted 
as compared to the linear methods. 

Providing a capability for the simulation of the nonlinear interaction of fluids and rigid body motion 
will be useful in several ways. First, the simulation could be used to validate the implementation of a 
control law derived using conventionally determined sensitivity coefficients. Although one must be wary 
of results obtained numerically, simplifications used to compute these force and moment derivatives 
will not be present in a Navier-Stokes simulation. Second, the coupled simulations could be used to 
develop a control law where nonlinear effects are important, using the computed aerodynamic forces and 
moments instead of tabulated empirical relationships. 5 Hence, computational design and prototyping 
of the aircraft control system offers the capability of reducing aircraft design cycle cost and enhancing 
safety as well as improving performance. 

The effort documented here begins to address the interaction of the disciplines of fluids, rigid 
body dynamics, and controls in nonlinear flight regimes. In order to assess the accuracy of these 
initial computations, a problem which could be compared against analytic and experimental results 
was chosen: the cavity store separation problem. Dix and Dobson's 6 recent experimental study of 
the separation of stores from cavity bays will be used as a basis for comparison. These wind-tunnel 
tests determined the trajectory of an uncontrolled missile from a rectangular cavity and will be used 
to validate the canard-fixed computation. 

Previous computational efforts have shown that the component problems of cavity flows 1 * 9 and 

uncontrolled store separation 10, 11 can be solved with reasonable accuracy using overset grid methods. 
It should be noted that inviscid solutions to the cavity store separation problem, with bodies in relative 
motion, have been obtained on tetrahedral meshes. 12 However, although these unstructured methods are 


2 



geometrically powerful. a means of consistently accounting for viscous effects is currently larking. Here, 
the combined problem of viscous cavity flow, rigid body dynamic's, and automatic' control techniques 
is addressed in a general manner using the overset mesh framework. 

The following sections discuss the approach used to solve the coupled system and tin? results ob- 
tained for several two- and three-dimensional cases. Comparisons of numerical results are made* against 
linearized or experimental results where available. 


Approach 

The procedure used to solve these coupled problems is shown schematically in Fig. 1. The process 


i 

I 

s 

35 

£ 



Fig. 1: Overall coupled system approach 


begins with the generation of the blocked mesh system while holding the relative position of the 
geometry in the initial position. After establishing the initial grid connectivity information, the solution 
of the flowfield can begin. Obviously, depending on the problem being addressed, this fixed grid solution 
may be steady or possess a bounded envelope. When this fixed-geometry flow solution is satisfactory, 
the motion of the body commences. 

The integrated aerodynamic loads acting on the body are provided to the G-DOF portion of the 
code along with the body state vector. Using this information, the body position and attitude is then 
integrated one time step, and kinematic constraints are applied. The controller also uses the body state 
to compute new effector settings, which are then applied to all relevant grids. Finally, since the grids 
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have changed relative positions. the intergrid communication is re-established. This process repeats 
for each time stop until the simulation is comploto. A detailed description of the compoiinit processes 
is given below. 


Flow Solver 

The fluid field was computed via the Reynolds-averaged Navier-Stokes (RANS) equations using the 
diagonal scheme of Pulliam and Chaussee 13 implemented in the overset, grid framework of Benek. Bull- 
ing, and Steger. 14 The equations were integrated through Euler implicit time marching and second-order 
spatial differencing with viscous wall conditions specified as no-slip, zero normal pressure gradient, and 
adiabatic. Information transfer across overset mesh boundaries was implemented using trilinear inter- 
polation of the dependent variable vector, Q = [p, pu, pv , pw, e ] T . The flow solver cost is 13/ns/cell/step 
on a Cray Y-MP. 


Turbulence Model 


In these computations, the slowly time-varying component of the flow is resolved, while rapid fluc- 
tuations are modelled. The algebraic turbulence model of Baldwin and Lomax, 15 as implemented by 
Buning, 16 is described below using a flow in the (x,y) plane. 

The wall-bounded flows use the original Baldwin- Lomax model, with the addition of a variable 
F max cutoff. The grids are chosen such that a unique wall distance is readily available. 

The cavity aperture spanning shear layer uses an eddy viscosity developed using F(y) = y |u;|, as 
suggested by Baldwin and Lomax 15 for wake regions. This results in 


Fwake 


oc 


C w 


wk~ 


ymax u dif 



max 



where specification of C w k is discussed below and the velocity difference is modified to be half the total 
velocity difference between the streams in the specified shear layer region 

u di f = J (u 2 + t/ 2 ) - J(v? + u 2 ) 

J V max V |u/| mai 


The free shear layer model is now given by 


Pt — pFC C pC w k 


V di£ 

|^| max 


(1) 


after dropping the Klebanoff intermittency function. Unmodified model constants K = 0.01G8 and 
C cp = 1.6 are used. 

The magnitude of the eddy viscosity in the free shear layer model can be altered with C w ^. Esti- 
mation of the proper value of C w k begins by using Gort.lers shear layer solution: 


u = 


U\ + U ‘2 
2 


1 + 


u 2 ~ 11 1 
112 + Wl 


erf(0 



where u\ and U 2 are the velocities of the slow and fast streams and £ is the similarity coordinate. The 
spreading parameter o is inversely related to the spreading rate, db/dx, where b is a measure of the 
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>hr;ir 1 ;iv<t width. I hr value of the spreading parameter when the velocity of one of tin- stream 
zero is rr, , . 

Gortler s solution ran be used to determine the maximum vorticity magnitude as follows: 
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Now. using PrandtFs mixing length assumption and scaling laws for jet boundaries, eddy viscosity 
can also be expressed as 


Pt oc pi 2 \uj\ 


= KopbAu (2) 

where Kc\ = 

4b(7~ 

Setting Eqs. 1 and 2 equal results in 

c w k = (ooKCcp^y 1 

and only cq remains to be specified. Empirical estimates of ctq range from 9.0 to 13 . 5. 9 For this series 
of cavity flow efforts ctq was set to 11.0, resulting in a value of C w k — 1 . 91 . 

Grid Generation and Communication 

Computation of the loads generated by cavity flows requires accurate representation of the geometry 
as well as the flowfield. Typically, a significant effort in grid generation is required before flow analysis 
can begin. Since matching zone faces are not required for the overset method used here, recent advances 
in algebraic 1 ' and hyperbolic 18 methods can be used. Hyperbolic grid generation, which gives good 
spacing and orthogonality control, was used for the walbbounded regions, while algebraic grids were 
used in shear flow regions. Advantages of this type of grid system include straightforward specification 
of the turbulent regions and allowance for independent refinement of each zone. This topology also 
permits the re-use of meshes for configuration studies. 

Exchange of flow information is accomplished using a domain connectivity function, with the donor- 
receiver relationship established at each time step using an efficient technique. 19 Although the cost of 
re-establishing intergrid communication is problem dependent, the computational expense is generally 
half of that used by the diagonalized flow solver. The initial location for the hole-cutters was specified 
using a graphical interface, 2 ® after which the movement of the grids and hole cutters were updated 
automatically. An example of the overset mesh topology used is shown in Fig. 2, which shows the 
configuration at an instant in the separation process. 


Kinematics and Rigid Body Dynamics 


Although details of the six-DOF dynamical equations can be found elsewhere, 21 briefly the rotational 
dynamics is described by Euler's equations of motion which aligns the xyz coordinates with the body 
principal axes at the center of gravity. For instance, for a rotation 9 about an axis A, Euler parameters 
can be specified as 
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Fig. 2: Coarsened grids after release 


These Euler parameters, integrated according to the rotational body dynamics, are updated and stored 
for each grid. 10, 19 Kinematic constraints can be imposed during the ejection process or for restricted- 
DOF simulations. In addition, the assumption of rigid-body dynamics eliminates the need to store the 
component grids for all time, since the Euler parameters may be used to compute grid attitude from 
the initial position. 

Rotating effectors were implemented by summing the relative commanded effector and body angular 
velocities, u. Integration of u) gives the proper effector attitude relative to the initial grid positions. 
The hinge line location is updated according to the attitude of the body. Storage of the hinge line and 
Euler parameters associated with each grid allows nesting of parent-child bodies to an arbitrary level 
without modification to the grid communication and support software. 

Pitch Attitude Control Law 

The fourth-order system shown in Fig. 3 was used as the system model for both the two- and 



Fig. 3: 2-D missile: Block diagram 

three-dimensional missile cases. In Fig. 3 the plant and servo can be expanded as 

0 _ ds + e 1 

S c as 2 + bs -I- c' r 

where the servo time constant was taken as r = X s. The plant coefficients were determined from 
the governing equation: they are composed of geometric and flow information, as well as the stability 
derivatives. The stability derivatives were determined from linearized supersonic airfoil theory in the 
two-dimensional cases and from direct computation in three-dimensions. Pole placement techniques 
and linearized system time response were used to determine the proportional. K a , and the rate gyro. 
AT, sensitivities. 
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Results and Discussion 


The method described above was applied in two-dimensions to a perturbed one degree-of-freedom 
case and a three degree-of-freedom cavity store separation process. The three-dimensional results are 
in progress. 


One Degree-of- Freedom Simulation 

In order to provide some measure of validation, a two-dimensional, 1-DOF simulation was imple- 
mented at a flight, altitide of 45,000 ft. Comparison of the linearized and the coupled nonlinear system 
response for a small perturbation allowed assessment of the basic methodology. 


Linearized Dynamics 

The equation of motion, M y = I yy 0, after neglecting the pitch-damping term, is 

My = qoo {[(C la Sd) c - (C la Sd) t }9 + (Ci a Sdf>) c } 

where the body was pinned at the location of the center of gravity for this analysis. 

Expressing the system in a state variable representation, with x = [6,(h0.6 c ] T , then the open loop 
equation can be expressed as 
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where 


x = Ax + Bu\ y = Cx ; u = Qx H- Rr 


The closed loop equation is expressed as 


x = Ax + B[Qx + Rr] = [A 4- bQ]x + K a br 
A + bQ — 


0 

0 
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1 
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0 
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-fKa -g{K a + \K r ) -gKr -f 
K a b — [0, 0. 0, (K a + K t )g} T 

The state can be computed using Euler explicit integration 


c n+1 = At 


(A + bQ + —)x n + K a br n 


; y 


n + 1 _ 


Cx 


n + 1 


The gain levels. K a = 1.2 and K r = 0.06 s , were computed using conventional root locus methods. 
These gains were input into a linearized one degree-of-freedom dynamics routine to verify implementa- 
tion of the control law. The result of the linearized analysis is shown in Fig. 4, which shows a slightly 
divergent envelope for the canard-fixed case, and damped behavior for the closed-loop system. Note 
that a canard deflection limiter was applied to represent stalling of the airfoil. 
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Fig. 4: Linearized 2-D missile: body orientation and canard deflection time histories for both uncontrolled 
and controlled cases 

Nonlinear Coupled Simulation 

The control law developed above was also implemented into the code which coupled the RANS 
equations with rigid-body dynamics. From Fig. 3, the control law can be expressed as 

6 C = — 7—7 \K a (6 c -9)- K r 9 + K,0 C } 

s + J L J 

and integrated using Euler explicit integration 

6 " +1 = /if+i { \ Ka(9c ~ sn) ~ Kr9n + Kt0c \ + *"} 

which can then be subjected to limit constraints. The result of the nonlinear solution shown in Fig. 5 
for both canard-fixed and controlled cases. After convergence to a steady-state solution, an accuracy- 
limited time step size of 4 Gfis was used, with a computational cost of five Cray Y-MP CPU hours. 
Grid communication was approximately one-third of the overall CPU time. Figure 5 shows similar 
behavior to the linearized cases, again with a modestly growing envelope for the canard-fixed case, and 
damped oscillations for the controlled case. This comparison provides a degree of validation of the 
implementation of the control law in the coupled code. 


2-D Resonating Cavity 

In order to establish confidence in the numerical method, a two-dimensional cavity computation 
was undertaken to demonstrate and validate self-induced cavity resonance. Details of this and three- 
dimensional cavity computations may be found elsewhere/ 

The test conditions, specified to match experiment, 22 were 

M 0 c = 0.9, Re L = 1.47 x 10 6 , L = 8 in. 

Poe — 0-40 kg/m 3 , p x = 2.9 x 10 4 N/m 2 
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Fig. 5: 2-D missile: Body orientation and canard deflection time histories and instantaneous C p contours 
for both uncontrolled and controlled cases at t = t\ 
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and a cavity length by depth. Lj D. of 2. Characteristic inflow and outflow conditions were >pecilied 
alony with a step size of Sf = 1.97/z.s. Power spectra were computed using 8192 ste]>s following rhe 
dissipation of the initial transients. 

Inspection of the computed cavity pressure history, shown in Fig. Ga. confirms the idealized feedback 



Fig. 6: 2-d cavity: (a) pressure history and (b) power spectra comparison 

mechanism identified from Rossiter’s experiments. Briefly, the cycle begins with the propagation of an 
acoustic wave from the aft wall of the cavity to the forward bulkhead. Wave reflection from the forward 
wall causes the shear layer to bow outwards, shedding vorticity. The deflected shear layer convects 
downstream and induces another cycle. This coupling of the acoustic and vortical fields is quantified 
by Rossiter's empirical model, given in Fig. 6b, which gives only feedback frequencies. 

In the frequency domain, comparison of Rossiter’s data to present results indicate agreement in 
frequency at the peak magnitudes, as shown in Fig. 6b. Magnitudes are higher for the present case by 
about 2 dB. which can be explained from dimensionality arguments. The solution was also found to 
be insensitive to second-order dissipation levels within the range 0.3 to 0.5. Figure 6b also shows the 
resonant modes predicted by Rossiter's equation, showing that K = 0.56 gives better prediction of the 
higher modes. Finally, the vertical knife edge schlieren images of Fig. 7 show the qualitative agreement 
between computed and observed 23 radiation patterns. 
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Experiment, Karamcheti Present 

Fig. 7: Comparison of experimental and numerical schlieren images with knife edge vertical 


2- D Controlled Store Separation 

Using the control law developed for the 1-DOF simulation above, albeit with an arbitrarily chosen 
9 C = 5°, a 3-DOF simulation of a store separating from a cavity was implemented. The separation 
began with the application of an 18000 N ejection force for 0.044 s, with corresponding acceleration 
of about 20#, during which the controller was off and no angular velocity was imparted to the store. 
The solution was initialized with the store fixed in carriage position, allowing damping of the starting 
numerical transients. Following convergence of the cavity acoustic envelope, the same accuracy-limited 
time step size of 4G [is was used. The time step size was chosen such that the streamwise Courant 
number was about unity in the shear layer. This restriction is equivalent to allowing an acoustic wave 
to propagate only one cell in a single step. Ten grids were used for this 1.3 x 10 5 point domain. The 
result, shown in Figs. 8 and 9, shows that the nose of controlled store remains pointed away from the 
parent body, while the canard-fixed store is pointed towards the parent 0.3 s after release. Since the 
controlled store is commanded to point away from the parent, the separation is faster for the compiled 
case than for the canard-fixed store. Inspection of the normal force history shows a component of about 
50 Hz. corresponding to the second stage of Rossiter's formula. 9 

3- D Missile Stability Derivatives 

In order to determine the proper feedback gains, the stability derivatives of the missile must be 
computed. For this 3-DOF simulation this includes C ma ,C m ^, and C mq . These parameters were 
computed from four cases: using a nominal 9 = 6 C = 0, a constant pitch attitude 9 = 0.04 rad . = 0, 

a constant deflection angle 9 = 0, 6 C = 0.04 rad, and a constant pitch rate 9 = 9 rad/ s. Canard motion 
was permitted by the yg in. nominal gap between the missile body and canard. This is a similar 
arrangement to the actual missile, albeit without the connector pin in the numerical model. The force 
and moment history was converged three orders of magnitude, approximately 2000 steps, oil this 18 
grid solution containing approximately 1.5 million points. Figure 10 shows the geometry used for this 
portion of the study, along with coefficient, of pressure, C p , contours. 


3-D Cavity Store Separation 

The geometry, shown in carriage position, can be seen in Fig. 11. The domain contains about 2.2 
million points distributed in 25 grids, with the missile grids being re-used from the previous stability 
derivative study. An example of the initialization of the cavity store separation problem is shown 
in Fig. 12. This store-fixed simulation will be run approximately five characteristic times, until the 
artificial starting transients have dissipated, after which ejection will begin. The ejection forces applied 
to the store will be such that the velocity at the end of the 8 in. piston stroke was 30 ft/s normal to 
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t=t 0 +79ms 




Fig. 8: 2-D controlled store separation: instantaneous Mach contours 
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Body orientation, 9 
Canard deflection, 6 


Body Attitude and 
Canard Deflection 


Time, t, s 


Fig. 9: 2-D controlled store separation: time histories of a) center of gravity, b) force and moment 
coefficients, and c) body attitude and canard deflection angles 



Fig. 10: 3-D missile: C p contours 
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the cavity with a pifch-dnwu rate of 1 rad/s. These* rates arc the nominal conditions used durum; wiud 
tunnel testing. 



Fig. 11: 3-D store separation: grids for missile in carriage 
Conclusions 

A pitch attitude control law was implemented in a coupled Navier-Stokes/rigid-body dynamics code. 
Comparison of nonlinear with linearized results showed reasonable comparison for the perturbed case 
with the controller active or off. Application of the control law to a two-dimensional, three-degree-of- 
freedom cavity store separation revealed improved trajectory characteristics. The generalized coding of 
aerodynamic effector kinematics in the coupled code will allow rapid implementation of existing control 
laws. Simulation of the coupled nonlinear aircraft trajectory can then be used to computationally 
prototype the control system. 


Note to the Reviewer 

Computation of a three-dimensional canard-fixed store separating from a cavity is in progress and 
will be compared to the experimental data of Dix and Dobson . 6 Comparison of the numerical and 
experimental results will provide an additional measure of validation. After completion of the uncon- 
trolled case, a simulation of the pitch-controlled case will commence. Determination of the feedback 
gains for this three-dimensional simulation will require the computation of the stability derivatives, 
which will also be accomplished via the Navier-Stokes equations. Assessment of the uncontrolled and 
controlled cases will show if improved cavity store separation characteristics can be achieved via canard 
effectors. 
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Fig. 12: 3-D store in carriage: instantaneous contours of pu on the symmetry plane and C }) contours on 
the wetted surfaces 
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